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ABSTRACT 


A complete  description  is  given  of  the  design,  implementation  and  use  of  a 
family  of  very  fast  and  efficient  large  scale  minimum  cost  primal  network  programs. 
Choice  of  data  structures  and  computational  testing  of  the  network  system  GNET  are 
discussed.  Important  extensions  are  explained  such  as  exploitation  of  special 
problem  structure,  element  generation  techniques,  post  optimality  analysis,  operation 
with  problem  generators  and  external  problem  files,  and  generalization  beyond  pure 
network  models. 


INTRODUCTION 


This  paper  reports  the  development  of  a large  scale  primal  network  code  that  is 
possibly  the  fastest  and  most  efficient  program  currently  available  for  solving 
capacitated  transshipment  problems.  The  capacitated  transshipment  problem  is  the 
most  general  of  the  minimum  cost  flow  models  which  include  the  capacitated  and 
uncapacitated  transportation  problems  and  the  personnel  assignment  problem.  These 
models  are  used  for  a large  number  of  diverse  applications  that  include  transporta- 
tion of  goods,  design  of  communications  and  pipeline  systems,  assignment  of  men  to 
jobs,  bid  evaluation  and  production  planning.  For  further  discussion  of  these 
applications  see  the  survey  articles  of  Bradley  [5],  Elmaghraby  [17]  and  Fulkerson 
[21]  and  the  textbooks  of  Busacher  and  Saaty  [8],  Charnes  and  Cooper  [9],  Dantzig 
[14]  and  Ford  and  Fulkerson  [19]. 

The  capacitated  transshipment  model  and  its  specializations  are  minimum  cost 

network  flow  problems.  The  goal  is  to  determine  how  (or  at  what  rate)  a good  should 

flow  through  the  arcs  of  a network  to  minimize  shipment  costs.  The  network  is  a 

directed  graph  defined  by  a set  of  nodes,  N,  and  a set  of  arcs.  A,  with  ordered 

pairs  of  nodes  (tail,  head)  as  elements  indexed  by  k.  For  each  arc  there  is  a 

shipping  cost  per  unit  flow,  c^^,  a minimum  allowable  flow  (or  lower  bound), 

and  a maximum  allowable  flow  (or  capacity)  , Each  node  is  either  a supply  node 

where  units  of  the  good  enter  the  network,  a demand  node  where  units  leave,  or  a 

transshipment  node.  The  problem  is  to  minimize  total  costs  with  flows,  x^,  that 

satisfy  the  associated  lower  bounds  and  capacities  and  preserve  the  conservation  of 

flow  at  each  node:  min  / 

keA 

s.t.  [ \ ~ I » i e N (1) 

keA  with  keA  with 

tail  i head  1 


u 


where  = { supply  if  1 is  a supply  node;  -demand  if  i is  a demand  node;  0 
otherwise  }. 

We  choose  this  notation  to  emphasize  that  data  will  be  stored  only  for  arcs  that 

are  present  in  the  network,  and  operations  are  defined  only  for  use  with  these  arcs. 

The  usual  textbook  notation  with,  for  instance,  I ^ x is  particularly  inappropriate 

i j 

since  for  large  practical  problems  it  is  never  true  that  all  node  pairs  are  connected 
by  an  arc.  I'urther,  our  notation  allows  multiple  arcs  (connecting  any  given  pair  of 
nodes)  which  are  common  in  practical  problems.  This  notation  is  also  consistent  with 
the  input  requirements  of  all  contemporary  large  scale  network  optimization  and  linear 
programming  codes. 

These  models  can  be  solved  as  linear  programming  problems  with  a constraint  for 
each  node  and  a variable  for  each  arc.  For  large  scale  problems,  contemporary 
commercial  linear  programming  codes  require  50-200  times  as  much  computer  time  and 
considerably  more  space  for  data  storage  than  special  purpose  network  flow  algorithms. 

The  transportation  model  was  originally  proposed  by  Hitchcock  [36]  1941  and 
Koopmans  [40]  1946.  Both  presented  computational  methods  that  would  now  be  called 
"primal  simplex."  Hitchcock  showed  that  an  optimal  solution  would  be  an  extreme 
point  solution  and  showed  how  to  iteratively  construct  better  extreme  point  solutions. 
He  noted  alternate  optimal  solutions  and  degeneracy  (all  in  6 1/4  pages).  Koopmans 
developed  simplex  multipliers,  or  "node  potentials"  and  the  optimality  criterion  and 
he  showed  that  an  extreme  point  is  equivalent  to  a tree.  Dantzig  [13]  1951  showed 
how  the  transportation  problem  could  be  solved  by  his  simplex  algorithm;  he  also 
developed  a special  variant  of  the  simplex  algorithm  for  the  transportation  problem. 
Orden  [48]  showed  that  these  results  can  be  extended  to  the  transshipment  problem. 

Network  models  are  widely  used  because  they  accurately  describe  a variety  of 
important  applications,  and  network  algorithms  can  efficiently  solve  large  problems 
(many  thousands  of  equations  and  variables) . Such  models  are  popular  because  they 


are  more  readily  accepted  by  nonanalysts  than  perhaps  any  other  class  of  opera- 


r 


tions  research  models.  Network  algorithms  can  economically  solve  problems  with  more 
variables  than  virtually  any  other  optimization  technique.  There  has  been  a surge  of 
interest  in  these  models  because  more  efficient  computer  programs  have  made  possible 
the  economic  solution  of  much  larger  problems. 

Although  many  papers  have  been  written  in  this  general  area,  and  significant 
computational  breakthroughs  have  been  reported,  there  has  not  previously  been  a single, 
unified  description  of  a complete  implementation,  nor  have  "new  generation"  computer 
programs  been  made  generally  available,  here  we  report  the  research  and  computational 
experiments  t^hich  have  produced  GNET,  an  extremely  efficient  yet  relatively  simple 
code.  An  important  objective  of  this  paper  is  to  make  these  new  approaches  easily 
accessible  to  a wide  audience  via  a clear  mathematical  exposition  and  a concrete 
example  of  a highly  efficient  FORTRAN  program.  Further,  the  availability  of  the 
computer  program  will  now  make  it  possible  for  any  investigator  to  reproduce  and 
extend  the  experimental  results. 

The  approach  we  have  used  in  this  research  is  to  study  the  data  structures  which 
seem  to  be  most  fundamental  in  the  sense  that  they  can  be  applied  to  many  types  of 
mathematical  programming  situations.  In  this  context,  we  view  the  major  advances 
over  the  last  thirty  years  in  efficient  solution  of  large  linear  programming  problems 
(for  example:  revised  simplex,  product  form  inverse,  bounded  variables,  generalized 

upper  bounding,  factorization,  sparse  matrix  methods,  etc.)  as  major  changes  in  the 
representation  of  the  data  accompanied  only  by  necessary  modifications  of  the  simplex 
procedure  to  accommodate  these  new  data  structures.  The  computational  breakthroughs 
in  primal  network  codes  are  also  due  to  improvements  in  data  representation  and  renew 
our  interest  in  the  subtle  relationships  between  the  algorithms  and  the  data  struc- 
tures employed  for  implementation. 

It  is  helpful  to  view  networks  as  an  important  special  case  of  large  scale  linear 
programming.  This  approach  is  crucial  and  overdue  because  there  has  been  very  little 
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success  in  extending  the  graph  theoretic  basis  tree  approach  to  more  general  mathe- 
matical programming  models.  For  instance,  the  graph  theoretic  proofs  of  the 
mathematical  correctness  of  generalizations  of  pure  network,  minimization  problems  are 
arduous.  Thus,  in  this  respect  we  consider  a purely  graph  theoretic  approach 
to  be  a dead  end.  We  have  taken  an  entirely  different  tack  in  developing  the  mathe- 
matical framework  of  a general  linear  programming  problem  and  specializing  it  to 
primal  network  models.  Therefore,  it  is  necessary  to  build  a sufficient  mathematical 
foundation  to  answer  the  question:  "What  is  the  capacitated  transshipment  problem  a 

specialization  of?",  rather  than,  "How  can  we  generalize  basis  trees?"  Thus  we  invite 
the  reader  to  view  this  paper  as  drawing  from  a general  theory  of  large  scale  mathe- 
matical programming  for  which  the  network  algorithms  are  concrete  realizations  of  a 
rich  algebraic  view  of  a problem  with  special  structure. 

Approaches  other  than  the  primal  simplex  have  been  proposed  including  out-of- 
kilter  (Fulkerson  [20]),  primal-dual  (Ford  and  Fulkerson  [18]),  dual  (Balas  and  Hammer 
[1]),  path  (Busacker  and  Gowen  [7]),  negative  cycle  (Klein  [38])  and  scaling  (Edmonds 
and  Karp  [16]).  Also,  special  algorithms  have  been  developed  for  the  assignment, 
shortest  route  and  maximum  flow  problems. 

The  contemporary  implementations  of  the  primal  simplex  network  algorithm  are 
based  on  compact  representation  of  the  basis,  determination  of  the  outgoing  arc  without 
trial  and  error  and  efficient  techniques  to  update  the  simplex  multipliers  at  each 
pivot.  Some  of  these  developments  were  suggested  in  the  1960's.  Glicksman,  L. 

Johnson  and  Eselson  [25]  described  a transportation  problem  with  few  sources  and  many 
sinks;  their  data  structure  may  be  viewed  as  storing  the  basis  as  an  arborescence  and 
using  this  structure  to  efficiently  find  the  outgoing  arc.  E.  Johnson  [37]  described 
a triple  label  scheme  that  represents  the  basis  as  an  arborescence  and  allows  the 
simplex  multipliers  to  be  efficiently  updated.  Johnson  describes  his  work  as  a modi- 
fication of  the  proposal  of  Scoins  [51] . 
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The  proposals  of  Glicksman,  L.  Johnson  and  Eselson  [25],  Scoins  [51]  and 
E.  Johnson  [37]  for  primal  algorithms  were  not  immediately  pursued;  the  most  widely- 
known  code  of  the  decade  was  an  out-of-kilter  implementation  by  Clasen  [12].  The 
contemporary  work  on  network  algorithms  was  begun  in  1970  by  Srinivasan  and  Thompson 
[54,  55],  Glover,  Karney  and  Klingman  [26]  and  Glover,  Karney,  Klingman  and  Napier 
[27].  This  work  was  a break  with  the  past  in  that: 

1.  Primal  algorithms  were  considered  despite  all  the  experiments  in  the  1950 's  and 
early  1960 's  that  showed  the  apparent  superiority  of  the  out-of-kilter  algorithm, 

2.  Contemporary  computer  science  tools  that  had  not  been  available  a decade  earlier 
were  used,  and 

3.  Computer  codes  were  developed  for  much  larger  problems. 

Later  and  independently,  McBride  [44]  and  Graves  and  McBride  [32]  specialized  their 
work  on  factorization  of  linear  programs  to  transshipment  problems.  Although  their 
development  was  quite  different,  the  network  specialization  of  their  data  structures 
is  similar  in  many  respects  to  data  structures  that  evolved  from  a graph  theoretic 
view  of  networks.  Mulvey  [45]  has  developed  an  efficient  large  scale  primal  code  at 
TRW.  Harris  [33]  has  developed  a primal  algorithm  for  transportation  problems  with 
many  sinks  and  few  sources.  Langley,  Kennlngton  and  Shetty  [42]  have  also  developed 
a primal  transshipment  code. 

A significant  ispect  of  contemporary  network  research  has  been  the  computational 
testing  of  different  algorithms  on  large  standard  test  problems.  One  major  topic  has 
bec.T  primal  algorithms  versus  out-of-kilter  algorithms.  Experiments  in  the  1950's  and 
early  1960 's  convinced  researchers  that  the  out-of-kilter  algorithm  was  superior, 
especially  for  transshipment  problems.  The  most  comprehensive  recent  comparison  has 
been  done  by  Glover,  Karney  and  Klingman  [26]  and  Barr,  Glover  and  Klingman  [2]  who 
compare  the  algorithms  on  a diverse  set  of  test  problems  [39].  Their  primal  code  was 
from  30%  (for  transshipment)  to  40%  (for  transportation)  faster  than  their  out-of-kilter 
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; code.  The  success  of  the  primal  algorithm  has  been  independently  verified  by  the 

experiments  of  others.  Most  researchers  now  believe  that  the  primal  algorithm  is 
, superior  to  others  including  the  out-of-kilter  algorithm  (an  exception  is  Hatch  [35]). 

Current  primal  implementations  are  faster,  require  less  storage,  are  more  suitable  when 
using  secondary  storage  devices  and  are  compatible  as  embedded  parts  of  more  general 

’ 

optimization  systems. 

Although  the  algebraic  development  of  our  computer  code  preceded  the  choice  of 
■ data  structures  and  the  algorithm,  we  can  also  establish  the  graph  theoretic  interpre- 

tation for  the  pure  network  problem.  Thus,  although  our  derivation  was  dissimilar 
from  historical  work  in  this  field,  we  are  able  to  show  the  relationship  'f  our  work 
to  that  of  others.  For  expository  purposes  we  will  draw  from  both  linear  algebra  and 
graph  theory,  using  pictures  and  terminology  consistent  with  past  literature;  the 
’>  detailed  algebraic  development  is  given  in  a companion  report. 

We  continue  now  with  a brief  algebraic  description  of  the  general  bounded 

i 

variable  simplex  algorithm  and  several  commonly  used  implementation  options.  The 
algebraic  specialization  of  the  simplex  method  for  pure  network  problems  is  presented. 
After  this  necessary  but  somewhat  mathematically  involved  section,  the  specific  design 
decisions  and  experiments  carried  out  with  GNET  are  described,  including  computational 
evidence  which  indicates  that  the  code  produced  by  this  approach  compares  very  favor- 
ably with  other  algorithms,  proprietary  or  otherwise.  Several  extensions  of  GNET  are 
presented,  including  codes  tailored  to  capacitated  and  uncapacitated  transportation 
problems,  and  other  variants  to  exploit  special  network  structure.  Postoptimal  and 
reoptimization  procedures  using  GNET  are  discussed.  Finally,  a review  of  the  litera- 
ture traces  the  original  contributions  found  to  be  of  fundamental  importance  in  this 
work. 

THE  PRIMAI,  SIMPLEX  ALGORITHM 

In  this  section  we  briefly  review  the  mathematical  underpinnings  of  the  bounded 
variable  revised  simplex  method.  In  order  to  maintain  a broad  scope,  the  presentation 
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intentionally  avoids  tangential  issues  of  specific  methods  and  implementations. 

Rather,  the  class  of  algebraic  algorithms  is  characterized  with  only  the  briefest 
indication  of  options  often  employed  for  actual  linear  programming  codes.  Small 
Roman  letters  will  denote  column  vectors,  and  a prime  indicates  transpose.  Large 
letters  denote  matrices;  superscripts  denote  columns  and  subscripts  denote  rows. 

Consider 

min  c'x  s.t.  Ax  = b and  0 ^ x ^ u ; (2) 

v/here  A is  viewed  for  the  present  as  a matrix  of  technological  coefficients. 

As  a practical  matter,  lower  bounds  on  the  variables  in  (1)  have  been  elimi- 
nated by  transformation. 

The  upper  bounds,  u,  are  most  readily  accommodated  implicitly.  Whenever  x^ 
reaches  its  upper  bound,  it  is  logically  replaced  by  Uj^  - Xj^;  column  k 

is  logically  treated  as  if  its  sign  has  changed  and  the  explicit  right  hand  side  is 
transformed  to  b - A u^.  If  a record  is  kept  of  each  variable  that  is  at  its  upper 
bound,  the  original  problem  solution  is  easily  recovered. 

By  construction  (possibly  including  introduction  of  unit  vectors  representing 
slack  and  artificial  variables)  A may  be  partitioned  A = (B,N),  with  B an  m x m 
matrix  of  linearly  independent  columns  which  is  called  a basis. 

Given  a feasible  basis  there  always  exists  a unique  x such  that 

BR  = b . 


In  terms  of  this  R there  is  always  a current  basic  solution  x°  = 
partitioning  c in  the  same  manner  as  A,  one  has 


(3) 

and,  upon 


= c-R 


(4) 


Any  generic  solution  x satisfying  the  constraints  can  be  rewritten  in  terms  of  the 
basic  solution: 


Ax  = (B,N)  = BXg  + NXj^  = b . 


(5) 


F 


Further,  since  B is  a basis,  there  exists  a unique  transformation  Z such  that 

BZ  = N . (6) 

Multiplying  (6)  by  and  subtracting  from  (3)  yields 

B(&-Zx^)  + Nxj^  = b , (7) 

and  subtracting  (7)  from  (5)  yields 

B(Xg-U-Zx^J)  = 0 . 

Since  the  columns  of  B are  linearly  independent,  x^  = - Zx^  and  the  general 

solution  becomes 

rx-zx^,-> 

(8) 


'X-Zx^- 

i ^ 

With  this  form  it  is  easy  to  compare  the  value  of  x to  any  current  solution 
x°  and  identify  an  improved  solution  when  one  exists.  The  value  of  the  generic 
solution  (use  (4)  and  (6))  is 

= c’x  + (c;-c^Z)Xj^ 

= + (Cj^-u'N)x^  ; (9) 

where  u (often  called  the  dual  solution  or  simplex  multipliers)  is  the  solution  of 

u'B  = c^  . (10) 

From  ( 9 ) it  is  clear  (since  x^^  ^ 0)  that  a necessary  condition  for  an  improved 
solution  is  that  there  exist  a column  of  N,  N , such  that 

c^  - u'N^  < 0 . (11) 

With  such  a column  consider  a specific  vector  obtained  from  (8)  by  taking  all 
components  of  x^^  zerc  except  Xj^, 


x — ( St"~Z  x^ , 0 , • • . , Xj^ , • • • ,0) 


(12) 
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As  a function  of 


^ 0,  this  vector  must  by  its  derivation  satisfy  the 


explicit  constraints.  Ax  = b ; feasibility  also  requires  satisfaction  of  the  bounds, 

0 ^ X ^ a.  For  components  > 0 this  requires  ^ 0>  or 

\ ^ V^ik  • ^13] 

For  components  z^,_  < 0 this  requires  x^  - ^ ^1’ 

(W/-==ik  • 

If  the  least  bound  on  Xj^  is  Uj^,  then  Xj^  stays  out  of  the  basis  but  goes  to 

its  upper  bound  and  (5)  and  (12)  yield  x=jc-Zu^  as  a new  basic  solution  with 

. .,k 

Bx  = b - N 

If  (13)  is  the  least  bound  on  Xj^,  taking  Xj^  = with  z^j^  >0  in  (12) 

d.  k. 

leads  to  the  exchange  of  B and  N in  the  partition  between  basic  and  nonbasic 
columns  and  the  new  basic  solution 

X = & - z , (x.  /z.,  ) , r i ; 

r r rk  1 ik 

and  X.  = x./z.,  • 

1 1 ik 


If  (14)  is  the  least  bound  on  x^^,  taking  Xj^  = with  z^^^  <0  in 

X k 

(12)  requires  the  basis  exchange  of  B and  N and  yields  x^  = x^  - z^j^(u.^-5i^) /-z^j^ , 

r i;  and  x.  = (u. -k. ) /-z ; as  a new  basis  with  x.  at  its  upper  bound. 

1 1 1 ik  1 

Assume  that  there  is  a current  basis  B,  a current  solution  k to  Bx  = b , and 

a current  solution  u of  u'B  = c' . A step  of  the  simplex  procedure  is  summarized: 

b 

k 

51.  Pfu.ce.0ut.  Select  a candidate  variable  to  enter  the  basis  with  (c^-u'N  ) < 0. 

k k 

52.  Pcotic  TcAt.  Find  the  greatest  bound  such  that  (with  BZ  = N ) : 

b)  for  z^^  > 0, 

c)  Xj^^  (u^-x^)Z-z^j^  for  z^j^  < 0. 


r 


If  the  minimum  ratio  is  determined  by  case  b)  or  c) , let  i be  a basis  variable  for 
which  the  mini  .um  is  achieved. 


S3.  P-i-VOt.  Update  the  primal  solution:  x = ^ - Xj^N  . If  Xj^  is  bounded  by  case 

a) ,  reflect  x.^  and  leave  the  basis  and  dual  sol'  -ion  unchanged.  Otherwise, 

i k 

change  th“  basis  by  exchanging  B and  N , for  case  c)  reflect  x^.  For  case 

b)  or  c)  find  the  new  dual  solution  to  u'B  = '. 

o 

In  executing  the  simplex  algorithm  a number  of  options  have  customarily  been 

k k 

employed  for  generating  the  solutions  of  the  linear  systems  BZ  = N and  u'B  = c'. 

B 

In  general  algorithms,  the  basis  inverse  Q = B ^ is  usually  used,  stored  and 
updated  in  some  form.  Further,  although  there  is  no  difficulty  in  deriving  a new 
algebraic  solution  to  (10),  u'B  = c^  , as  a practical  matter  u may  be  directly 
achieved  from  u by  simple  transformation. 


. th 


(15) 


Proposition : u = u + AQ^  where  is  the  i — row  of  the  inverse  of  B. 

k k k k Ic  If 

The  new  basic  column  N determines  A as  uN  = uN  + A(Q^N  ) = uN  + A(Q^BZ  ) 

k k i 

uN  + ~ ^k  ’ that  A = (Cj^-uN  . Exclusive  of  the  outgoing  column  B , 

u'E*^  = ub”^  + A(Q.B  ) = uB=c,  r ^ i.  □ 

1 r r r 

k 

The  (pivotal)  update  of  Q after  exchange  of  B and  N is  easily  derived 
(e.g.,  [30]).  The  most  elementary  and  explicit  procedure  is  to  carry  and  update  (by 
pivoting)  a complete  tableau 

r QN  Qb^ 

[c^-u’N 


The  revised  simplex  procedure  generates  these  elements  as  needed  by  access  to  columns 
of  N,  c and  b and  use  of  Q.  Most  full  scale  systems  employ  an  additional 
refinement  by  expressing  Q as  the  product  of  elementary  "eta"  column  vectors,  each 
representing  the  pivotal  transformations  generating  Q from  an  initial  basis. 
Frequently  the  history  of  'eta'  columns  grows  too  long  for  reasonably  efficient 
generation  of  Q,  or  numerical  error  is  propagated  and  detected,  forcing  a reln- 
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version  with  'eta'  column  representatives  from  only  the  current  basis. 

Other  systems  support  the  solution  procedure  by  using  combinations  of  features 
such  as  an  "LU"  decomposition  of  B [3],  a Cholesky  decomposition  of  BB'  = LL' 

[23,  49,  50],  or  list  representation  of  nonzero  elements  of  problem  components  and 
coefficient  representation  by  pointers  to  a table  of  real  values.  Hybrid  schemes 
factor  B into  partitions  with  special  structure:  Generalized  Upper  Bounding  (GUB) 

identifies  an  inherent  identity  matrix  for  some  rows  of  B [15];  a partial 
triangulation  of  B with  an  inverse  for  remaining  columns  and  rows  can  be  used  [31, 
32].  Whether  systems  solve  (6)  and  (10)  by  triangular  substitution,  inverse  trans- 
formation, or  some  combination,  all  are  algebraically  equivalent  simplex  implemen- 
tations differing  only  in  the  structures  chosen  to  support  computation  for  the  class 
of  problems  at  hand. 

PRIMAL  ^■ETW0RK  SPECIALIZATION 

A specialization  of  the  simplex  algorithm  to  the  transportation  problem  was 
developed  by  Dantzig  [13]  in  1951.  It  is  not  surprising  that  the  transportation 
algorithm  was  developed  immediately  after  the  simplex  algorithm,  because  the  works  of 
Hitchcock  [36]  1941  and  Koopmans  [40]  1946  on  the  transportation  problem  contain  many 
concepts  that  presage  the  simplex  algorithm.  The  interaction  of  general  linear 
programming  algorithms  and  transshipment  algorithms  has  a long  history  that  has 
enriched  the  study  of  both. 

Here  we  establish  explicitly  the  relationship  between  the  general  primal  simplex 
algorithm  and  the  modem  implementations  of  the  transshipment  algorithms.  Our  goal  is 
to  understand  the  algebraic  foundations  of  the  modem  transshipment  implementations. 
Also  (and  perhaps  more  importantly)  we  lay  the  groundwork  for  the  next  stage  in  the 
interplay  of  these  models:  the  incorporation  into  the  next  generation  of  general 

linear  programming  computer  systems  of  the  important  ideas  that  have  made  possible  the 
breakthroughs  for  transshipment  problems. 

The  fundamental  fact  that  permits  design  of  efficient  primal  transshipment 
algorithms  is  the  well-known  result  that  any  transshipment  basis  can  be  put  in  (upper) 
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r 

triangular  form  by  a simple  permutation  triangulation.  This  inherent  triangularity 
can  be  exploited  by  network  specializations  of  the  simplex  method  by  directly  solving 
(6)  by  back  substitution  and  (10)  by  forward  substitution.  Also,  the  triangulated  » 
basis  simplex  algorithms  lead  to  mucn  more  efficient  network  solutions  due  to  other 
fortunate  simplifications.  The  most  remarkable  of  these  is  that  the  solution  update 
of  step  S3  can  be  accompanied  (in  fact  aided)  by  a very  simple  and  efficient  dynamic 
re  triangulation  of  each  new  basis. 

An  initial  transshipment  basis  with  full  row  rank  can  always  be  constructed  by 
introducing  for  each  row  in  (1)  a unit  vector  with  sign  matching  that  of  the  right 
band  side  (negative  for  demand  nodes).  With  the  addition  of  these  artificial  vectors, 

A has  full  row  rank  and  each  column  of  A has  a single  nonzero  element  1,  a single 
nonzero  element  -1,  or  two  nonzero  elements  (a  1 and  a -1). 

Theorem:  Any  basis  B extracted  from  A for  a transshipment  problem  can  be 
triangulated  by  rearranging  rows  and  rearranging  columns. 

Proof:  Let  B have  m rows.  Locate  a column  with  a single  nonzero  entry.  Exchange 

rows  and  columns  so  that  the  nonzero  entry  is  the  first  diagonal  element.  For  the 
k.-^  step  of  the  construction  rows  and  columns  have  already  been  rearranged  so  that 
columns  l,2,...,k-l  have  only  zeros  below  the  diagonal.  Select  a column  with  a 
single  nonzero  entry  in  rows  k,...,m.  Exchange  rows  and  columns  so  that  the  nonzero 
elenx'nt  is  the  k—  diagonal  element.  This  construction  is  well  defined  only  if 
there  is  a column  with  a single  nonzero  element  in  rows  k,...,m.  There  must  be  such 
a column  for  otherwise  each  column  would  have  no  nonzero  elements  or  exactly  one  +1 
and  one  -1  in  rows  k,...,m  and  the  sum  of  rows  k,...,m  would  be  the  zero  row  which 
would  contradict  the  assumption  that  B is  a basis.  □ 

A graph  can  be  defined  that  represents  the  transshipment  basis.  Let  I = {ij^,l2» 
...,i^}  be  a row  ordering  corresponding  to  a triangulated  B.  Associate  with  each 
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node  the  row  number  off  diagonal  element  in  the  k — column  of  the 

triangulated  basis;  if  there  is  no  nonzero  offdiagonal  element,  set  P(ij^)  to  m + 1. 
Define  a graph  with  nodes  l,2,...,m+l  and  let  (ij^,p(ij^))  k = l,2,...,m  be  directed 
arcs  from  i.  to  p(l,  ).  Since  each  node  i,  is  connected  to  a node  p(i,  ) = i 
with  h < k or  to  node  m + 1,  there  is  a directed  path  from  each  node  to  node  m + 1. 
The  graph  is  called  a fioote.d  CUibo^e^cmce  [4]  with  node  m + 1 the  fwot.  Ignoring 
the  orientation  of  the  arcs,  since  the  graph  is  connected  and  has  m + 1 nodes  and 
m arcs  it  is  a tree  (it  can  be  shown  [4]  that  this  definition  is  equivalent  to  the 
usual  definition  cf  a tree  as  a connected  graph  with  no  cycles).  In  the  computer 
science  literature  (e.g.,  [41])  the  term  tree  is  often  used  instead  of  rooted 
arborescence . Figure  1 is  an  example  of  a transshipment  problem  with  a basic  feasible 
solution  specified  and  Figure  2 has  the  basis  and  the  associated  arborescence.  Our 
pictorial  representation  with  the  "root"  at  the  top  rather  than  at  the  bottom  is 
fairly  standard. 

For  node  i,  p(i)  is  called  the  pfizdndUiiOH  of  i and  the  rooted  arborescence 
is  called  the  pK'.dtdUiiOfi  QKa.pk.  The  predecessor  graph  is  closely  related  to  but  not 
identical  with  tne  classical  result  of  Koopmans  [40]  that  the  arcs  of  a transshipment 
basis  form  a tree  over  the  nodes  of  the  problem.  The  classical  tree  preserves  the 
orientation  of  arcs  in  the  original  network  and  does  not  include  node  m + 1. 

The  predecessor  function  p(  ) is  a well  known  compact  way  to  represent  trees 
and  rooted  arborescences  and  has  been  used  in  network  algorithms  for  at  least  15  years. 
The  predecessor  graph  has  often  been  used  interchangeably  with  the  classical  basis 
tree.  It  is  important  to  distinguish  between  them  because  the  predecessor  graph  is  a 
data  structure  that  supports  the  computation  of  the  algorithm  and  the  orientation  of 
the  arcs  indicates  the  unique  direction  to  the  root  rather  than  a direction  in  the 
network.  Furthermore,  the  predecessor  graph  can  be  extended  to  triangular  bases  that 
have  no  underlying  network. 


i 
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Figure  1.  A Single  Commodity  Transshipment  Problem 
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Transshipment  Basis  Matrix  A Preorder  Triangulation 


P(ij^)  = [834882  1 13  2 326] 

The  Basis  Predecessor  Function 


^ Root 


Figure  2,  A Transshipment  Basis  (for  Figure  1) . 
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The  predecessor  function  p(  ) can  be  iterated  to  construct  the  unique  path 
from  any  node  to  the  root;  this  path  is  called  a backpcUh.  The  hmtdixitt  4tUCC.U-i>0fi£ 
of  a node,  if  any,  are  the  first  nodes  encountered  on  all  paths  except  the  backpath 
to  the  root  and  all  the  nodes  on  these  paths  are  called  the  ■l>acceJ>60fUt . Another 
characterization  of  the  successors  of  node  i is  that  they  are  all  nodes  whose 
backpath  to  the  root  includes  node  i.  A tree  can  also  be  represented  by  the 
immediate  successors  of  each  node;  however,  since  the  class  of  trees  that  arise  in 
network  problems  is  called  m-ary  with  0 to  ra  immediate  successors  for  each  node, 
this  is  more  difficult  to  maintain  dynamically  than  the  predecessor  function  which 
always  has  a single  unique  value  for  each  node  except  the  root. 

In  general,  there  are  many  different  triangulations  for  any  given  transshipment 
basis.  (Note  that  at  each  step  in  the  construction  there  may  be  several  choices  for 
the  next  column.)  However,  all  such  triangulations  yield  the  same  predecessor  func- 
tion and  graph  (where  the  ordering  of  successors  right  to  left  for  any  node  is 
immaterial).  Thus,  the  predecessor  graph  does  not  completely  represent  a triangula- 
tion without  additional  information,  namely  an  ordering  of  the  rows.  A mathematical 
development  of  these  properties  and  their  implications  is  given  in  the  companion 
paper. 

The  relationship  between  the  algebraic  view  of  the  simplex  algorithm  and  the 
graph  theoretic  view  of  much  of  the  network  literature  can  be  shown  by  describing  the 
operations  of  the  simplex  algorithm  and  the  triangulation  in  terms  of  the  predecessor 
graph.  A graph  theoretic  proof  of  the  triangulation  theorem  identifies  a node  with  each 
row  and  an  arc  with  each  column  with  two  nonzero  elements.  Also  included  is  the  root 
node;  columns  with  a single  -1  (or  1)  are  represented  as  an  arc  to  (from)  the  root 
node.  A trian^ulation  and  the  predecessor  graph  are  constructed  by  first  selecting 
the  root  node.  Select  any  arc  with  one  end  the  root  and  add  the  new  node  and  arc 
orienting  the  arc  toward  the  root.  For  the  k—  step  of  the  construction  k nodes 
and  k - 1 arcs  are  already  in  the  graph.  Select  an  arc  with  exactly  one  node 
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already  in  the  graph,  add  the  new  node  and  the  arc  to  the  graph  orienting  the  arc 
toward  the  old  node.  The  resulting  triangulation  is  characterized  by  the  order  in 
which  the  nodes  and  arcs  are  introduced. 


For  the  network  linear  program  each  node  is  identified  with  an  equality 
constraint  and  each  arc  k from  node  i to  node  j is  identified  with  a variable 
with  column  N having  a 1 in  row  i,  a -1  in  row  j and  O's  elsewhere.  For  the 
discussion  below  it  is  assumed  that  the  basis  B is  in  triangular  form  and  to 
simplify  notation  it  is  momentarily  assumed  that  row  i of  B corresponds  to  node  i 
for  all  i e W.  (Renumber  the  nodes  if  necessary  so  that  I = {1 ,2 , . . . ,m} . ) It  is 
further  assumed  that  all  the  diagonal  elements  of  B are  1;  if  B initially  has  a 
-1  on  the  diagonal,  the  reflection  Uj^  - Xj^  transforms  it  to  a 1.  This  may  be 
viewed  as  transforming  the  basis  variables  to  make  the  orientation  of  arcs  in  the 
predecessor  graph  the  same  as  the  orientation  in  the  Koopmans  basis  tree. 

Given  the  vector  of  simplex  multipliers  u,  the  priceout  formula  (11)  for  non- 
basic  arcs  (step  SI)  simplifies  to  c,  - u.  + u..  Thus  the  priceout  involves  only 
addition  operations. 

For  the  determination  of  the  arc  to  leave  the  basis  in  step  S2  the  system  of 
Ic  k.  k 

equations  BZ  = N must  be  solved  for  Z . This  calculation  is  described  by 

showing  how  to  solve  BQ'^  = e^  for  Q'^ . (Q^  is  thus  the  column  of  the 

inverse  of  B.)  Since  B is  triangular,  can  be  obtained  by  simple  backward 

solution:  the  m,m-l, . . . , j+1  elements  of  are  seen  to  be  0,  the  j—  element 

is  1.  Setting  the  j—  element  equal  to  0 puts  a 1 in  the  modified  righthand  side 

in  the  row  corresponding  to  the  off  diagonal  -1  (if  any)  in  the  j—  column  of  B, 

i til 

this  row  is  p(j).  As  before,  elements  j+1, . . . ,p(j)-l  of  are  0 and  the  p(j)— 

element  is  1.  This  continues  putting  1 in  the  p(p(j))  ,p(p(p(j)))  , etc.  elements  of 
until  a column  with  no  nonzero  offdiagonal  is  encoiontered.  Elements  back  to 
elements  ...,2,1  are  then  set  to  0. 
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Thus,  in  terms  of  the  predecessor  graph,  all  arcs  traversed  on  the  backpath  from 


node  j to  the  root  have  an  element  1 in  ; all  other  elements  are  0.  Therefore,  Q"' 
can  be  generated  directly  from  the  precedence  function  p(  ),  which  (with  I)  gives 


the  substitution  rules  for  the  back  solution. 


k k i i 

The  calculation  of  Z follows  immediately  since  Z = Q - Q . The  element  in 


Ic  i i 

a row  of  Z is  1 for  all  rows  with  a 1 in  Q alone,  -1  for  all  rows  with  1 in  Q 


alone  and  0 for  all  rows  with  0 in  both  or  1 in  both.  Since  the  nonzero  elements  of 
,k 


Z are  1 or  -1,  the  calculations  in  steps  S2b  and  S2c  also  involve  only  addition  and 
subtraction.  Further,  the  calculation  is  usually  reduced  enormously  by  the  extinguish- 
ment of  the  elements  common  to  both  i and  j backpaths. 


As  an  example,  consider  the  triangulated  basis  in  Figure  2 (noting  that  the  nodes 


k 9 9 

have  not  been  renumbered)  with  column  N associated  with  arc  (1,9),  BQ  = e , and 


BQ^  = e^. 


I = {8  1 7 5 4 3 2 11  9 6 12  10} 

p(  ) = {8  3 4 8 8 2 1 13  2 3 2 6} 

= [110000000000]' 
= [1  0 0 0 1 1 1 0 1 0 0 0]’ 

Z*^  = [0  1 0 0 -1  -1  -1  0 -1  0 0 0]' 


In  step  S3,  the  primal  solution  is  updated  using  Z . Also,  as  shown  by  (15), 


the  simplex  multipliers,  u,  can  be  updated  rather  than  calculated  at  each  step. 


th 


The  algebraic  characterization  of  the  update  applies  the  i — row  of  the  inverse  of 


,th 


B,  Q^,  where  the  outgoing  arc  is  the  SL — column  of  B.  The  characterization  of 


above  showed  that  an  element  of  is  nonzero  and  equal  to  1 only  if  the 


corresponding  basic  arc  is  traversed  on  the  backpath  from  node  J to  the  root.  It 


, th 


follows  that  is  all  O's  except  for  1 in  the  H — element  and  a 1 for  any  element 


j such  that  the  outgoing  arc  is  on  the  backpath  from  j to  the  root.  In  terms  of 


the  predecessor  graph,  is  1 for  node  S,  and  each  of  its  successors  and  0 else- 
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where.  Since  the  nonzero  elements  of  Z are  1 and  -1,  A in  (15)  is  either  plus  j | 

or  minus  - u^  + u where  k is  the  incoming  arc.  Since  is  all  O's  and  I's,  i 

j 

the  update  of  u is  accomplished  by  adding  A to  u and  to  the  u's  of  (only)  the  j 

K/  j 

successors  of  node  1.  i 

As  will  be  shown  in  the  next  section,  the  simple  additive  updates  of  the  primal  I 

solution  and  dual  multipliers  in  step  S3  are  actually  accomplished  in  a single  inte-  | 

j 

grated  process.  This  ("pivot")  process  simultaneously  performs  the  updates  while  | 

retriangulating  the  new  basis.  Fortunately,  the  retriangulation  is  also  inexpensive  | 

to  perform  on  a transshipment  basis  with  a single  new  column  violating  triangularity. 

IMPLEMENTATION 


The  design  of  large  scale  programming  codes  necessarily  involves  many  significant 
decisions  which  have  major  impact.  The  following  fundamental  principles  were  used  to 
resolve  design  questions  in  the  development  of  the  code  reported  here. 

1.  The  code  is  designed  for  large  scale  problems.  Even  though  experimental 
testing  will  be  confined  by  economic  considerations  to  problems  with  some 
arbitrary  maximum  size  (say,  10,000  equations),  the  design  decisions  should 
lead  to  a code  with  superior  large  scale  performance. 

2.  The  goal  is  a code  for  the  most  general  capacitated  transshipment  problem. 

While  problems  with  specializations  (e.g.,  uncapacitated,  transportation, 
assignment)  must  be  solved,  the  basic  code  will  not  be  tailored  to  these 
special  features.  In  addition,  no  special  numbering  of  nodes,  extensive 
preprocessing,  or  other  design  specificity  will  be  required  that  will  limit 
the  capability  of  the  code.  Efficient  solution  of  problems  should  not 
require  detailed  advance  knowledge  of  problem  structure  (for  instance,  a 
feasible  initial  solution).  Problems  with  multiple  arcs  will  be  accommo- 
dated. 
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3.  For  practical  problems  the  number  of  arcs  is,  in  general,  much  greater  than  the 
number  of  nodes,  m.  However,  the  problems  are  usually  sparse  in  the  sense 


that  the  number  of  arcs  seldom  approaches  the  maximum  number  m(m-l)  of  oriented 
node  pairs.  Thus,  it  is  significantly  more  expensive  to  store  information 
associated  with  each  arc  than  it  is  for  each  node  and  prohibitively  expensive 
to  store  a node-arc  incidence  matrix.  Practical  general  minimum  cost  network 
flow  problems  are  always  heavily  capacitated. 

4.  It  is  important  to  produce  a code  that  is  machine  independent  as  well  as 
efficient.  For  example,  machine  specific  features  such  as  assembly 
language,  use  of  particular  offline  mass  storage  devices,  storage  of  data 
using  bit  string  logical  vectors,  use  of  other  architectural  curiosities 

on  particular  machines  or  nonstandard  language  features  are  all  to  be  avoiced. 
The  language  used  is  basic  FORTRAN. 

5.  Where  feasible,  speed  of  execution  will  be  given  preference  over  economy  of 
bpcRC  for  data  storage. 

6.  Since  the  program  will  be  used  for  comparisons  of  various  data  structures  on 
a wide  variety  of  network  problems,  it  must  be  equipped  with  effective 
external  tuning  parameters.  While  some  tuning  is  possible  for  the  pivot 
mechanism,  the  pricing  scheme  invites  especially  close  scrutiny  for  tuning 
purposes . 

Once  the  design  of  an  efficient  network  code  is  chosen,  consideration  will  also 
be  given  to  additions  of  other  advanced  features  such  as  "in-core/out-of-core" 
operation,  implicit  arc  generators,  crashed  bases,  nonlinear  costs,  post  optimal 
analysis,  and  so  forth.  These  extensions  will  not  be  allowed  to  interfere  with  the 
basic  design  goals,  but  they  should  not  be  precluded  by  the  basic  design  decisions. 
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The  description  of  GNET  begins  with  discussion  of  arc  and  node  storage 
representation.  Next,  the  ratio  test  is  described.  The  pivot  is  explained  function- 
ally (with  a more  extensive  algebraic  derivation  left  for  the  companion  paper). 

Finally,  the  pricing  mechanism  is  examined.  This  order  is  chosen  (steps  S2  and  S3 
followed  by  SI)  to  faithfully  report  the  implementation  historically  and  to  lead 
smoothly  to  the  computational  performance  tests.  Hereafter,  notation  with  upper-case 
Roman  letters  indicates  a program  variable  and  addition  of  parentheses  denotes  an 
array.  For  instance,  the  predecessor  OAAoy  is  referred  to  as  P(  ). 

Since  there  will  be  many  arcs,  it  is  critical  to  minimize  the  stored  data 
describing  each  arc.  A typical  input  format  (for  example,  SHARE)  for  each  arc  is: 
tail,  head,  cost,  lower  bound  and  upper  bound  on  flow.  The  lower  bounds  are  removed 
by  transformation.  If  the  arcs  are  sorted  so  that  all  arcs  with  the  same  head  are 
stored  in  contiguous  space,  the  list  of  heads  can  be  replaced  by  a node-length  array 
whose  j—  element  is  the  location  of  the  first  arc  with  head  j.  Since  network 
models  have  many  more  arcs  than  nodes,  this  reduces  the  storage  requirements  for  the 
algorithm.  Thus,  the  network  is  stored  as  three  arc- length  arrays:  the  tails  T(  ), 

the  costs  C(  ) and  the  upper  bounds  (capacities)  CP(  );  also,  one  node-length  array 
is  used  , the  head  entries  H(  ) into  T(  ) . Positive  capacities  are  required  after 
transformation  of  lower  bounds  for  all  arcs — uncapacitated  arcs  have  capacity  set  to 
some  value  greater  than  the  total  supply.  Arcs  out  of  the  basis  at  their  upper 
bound  are  marked  with  a sign  bit  on  the  capacity  (-CP(  )). 

It  is  natural  to  associate  with  each  node  i (except  the  root)  the  unique  basic 
arc  that  connects  i to  its  predecessor.  It  is  convenient  to  have  the  basic  arcs 
oriented  the  same  as  in  the  predecessor  graph,  that  is,  from  a node  to  its  predecessor. 
This  is  accomplished  by  reflecting  arcs  as  necessary.  The  predecessor  array  P(  ) is 
marked  with  a minus  sign  for  all  arcs  that  have  not  been  reflected.  (Subsequently 
when  using  P(  ) we  will  assume  for  simplicity  that  it  is  positive.)  The  flow  on  arc 
(i,P(i))  is  stored  in  X(i). 
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In  step  S2  of  the  simplex  algorithm,  it  is  necessary  to  compute  capacity  minus 
flow  for  basic  arcs  with  <0.  It  is  convenient  to  use  a node-length  array  to 

speed  up  this  calculation.  One  approach  is  to  store  for  each  basic  arc  a pointer  to 
its  capacity  in  the  CP(  ) array.  Another  technique  is  to  store  the  capacity  rather 
than  the  pointer.  A third  method  is  used  in  GNET.  The  capacity  minus  flow  is  stored 


in  a node-length  array  CPX(  ) . 

The  simplex  multipliers  are  stored  in  a node-length  array  U(  ).  Figure  3 shows 


these  arrays  for  the  basis  in  Figure  2. 

After  the  incoming  arc  has  been  chosen  in  step  SI,  the  outgoing  arc  is  determined 

in  step  S2  and  the  data  structures  are  updated  in  step  S3.  Let  k be  the  incoming 

arc  going  from  node  i to  node  j.  The  possible  outgoing  arcs  correspond  to  the 

Ic  k 

nonzero  entries  in  Z . We  have  already  seen  that  Z is  zero  for  all  nodes  common 
to  both  i and  j backpaths.  Let  the  join  be  the  first  node  on  the  backpath  from 
i to  the  root  that  is  on  the  backpath  from  j to  the  root  (the  join  is  the  extin- 
guishment  point  for  Z ) . The  possible  outgoing  arcs  are  the  arcs  on  the  backpath 
from  i to  the  join  and  the  arcs  on  the  backpath  from  j to  the  join. 

It  is  critical  to  identify  the  backpaths  from  i and  j to  the  join  efficiently. 
The  trial  and  error  of  the  classical  "stepping  stone"  methods  of  most  textbooks  will 
clearly  not  suffice  for  any  but  trivial  problems.  We  discuss  four  methods  to  identify 
the  backpaths  from  i and  j to  the  join  node.  Note  that  only  predecessor  infor- 
mation is  available  to  the  program  at  each  node  and  that  our  data  structure  has  no 
global  view  of  the  arborescence  as  does  the  reader  of  Figure  2. 

The  most  naive  method  is  to  mark  in  some  way  all  the  nodes  on  the  backpath  from 
i to  the  root.  Then,  the  first  marked  node  encountered  on  the  backpath  from  j to 
the  root  is  the  join.  This  method  is  satisfactory  for  smaller  problems,  but  for 
larger  problems  it  is  more  efficient  to  avoid  the  iteration  along  one  complete  back- 
path  all  the  way  to  the  root  by  keeping  additional  information  about  the  tree.  Also, 
in  this  way  it  will  not  be  necessary  to  unmark  the  marked  nodes  before  performing  the 
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Figure  3.  GNET  Arrays  (for  Basis  in  Figure  2) 


”^Two  candidate  queue  arrays  also  used  by  GNET  are  omitted  for  clarity.  Actual 
simplex  multipliers  in  GNET  would  all  be  exactly  1,188  units  smaller  than 
shown.  This  difference  is  the  high  cost  associated  with  artificial  demand 
arc  13->8,  initially  computed  as  BM  = Nodes  x maximum | cost | . 
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next,  simplex  pivot. 

One  efficient  method  is  to  store  for  each  node  the  number  of  nodes  on  the 
backpath  to  the  root,  call  the  de.ptk  of  the  node,  D(  ) , in  the  tree.  With  depth  » 

information  available  it  is  possible  to  iterate  the  backpaths  synchronously  from  i 

. < 

and  j to  identify  the  join  without  iterating  either  backpath  past  the  join.  Depth 

i 

is  used  to  indicate  which  backpath  node  is  deeper  in  the  tree  and  should  be  iterated. 

When  both  backpath  nodes  have  matching  depths,  the  nodes  are  compared  for  equality. 

A match  indicates  the  join,  and  a mismatch  indicates  that  both  backpaths  should  be 
iterated  for  another  comparison. 

Another  efficient  method  similar  to  the  depth  approach  is  to  store  for  each  node 
the  number  of  nodes  in  its  subtree,  called  the  numbe/L  , NS(  ).  Starting 

with  nodes  i and  j,  the  backpath  node  with  strictly  fewer  successors  is  iterated. 

When  both  backpath  nodes  have  the  same  number  of  successors,  a match  indicates  the 
join  and  a mismatch  forces  iteration  of  both  backpaths. 

The  fourth  method  will  be  discussed  below. 

The  latter  three  methods  for  locating  the  join  look  only  at  arcs  that  are  on  the 
backpaths  to  the  join,  thus  it  is  possible  to  determine  the  outgoing  arc  while 
searching  for  the  join.  As  noted  above,  all  the  arcs  on  the  backpath  from  1 to  the 
join  are  the  +1  elements  of  Z and  all  the  arcs  on  the  backpath  from  j to  the 

join  are  the  -1  elements.  The  ratio  test  step  S2  is  then  simply 

1CP(k)  the  capacity  of  the  incoming  arc 

X(  ) for  arcs  on  the  backpath  from  i to  the  join 

CPX(  ) for  arcs  on  the  backpath  from  j to  the  join. 

The  computational  simplicity  of  this  ratio  test  is  the  rationale  for  the  reflection 
of  basic  arcs  and  the  adoption  of  CPX(  ). 

If  the  incoming  arc  is  out  of  the  basis  at  its  capacity,  then  step  S2  may  be  . 

viewed  as  increasing  flow  in  a fictitious  arc  from  node  j to  node  i with  the  same 
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capacity. 

A major  proportion  of  the  work  of  each  simplex  step  is  S3,  the  pivot.  In  this 
step,  the  entering  and  leaving  arcs  are  exchanged,  the  flows,  X(  ) and  CPX(  ), 
are  updated,  the  simplex  multipliers,  U(  ) , are  changed,  and  the  simplex  data  arrays 
P(  ) and  D(  ) (using  the  dcnth  mechanism  for  example)  are  modified. 

Within  this  step  the  coordination  and  sequencing  of  operations  are  critical  to 
the  efficiency  of  the  network  algorithm  because  the  manipulation  of  many  nodes  and 
heavy  use  of  the  simplex  data  structure  are  involved.  If  properly  done,  this  step  is 
the  elegant  central  part  of  the  code  that  can  be  executed  by  a computer  quickly; 
however,  the  explanation  will  be  somewhat  intricate — this  part  of  the  algorithm  is 
considered  by  many  to  be  the  "secret"  part  of  the  proprietary  codes. 

To  illustrate  a typical  pivot.  Figure  2 shows  the  entering  arc  (i,j),  join  and 

leaving  arc  (c,d) . Call  the  backpath  from  j to  c the  p-ivot  itzm.  Figure  4 in- 
cludes the  predecessor  graph  after  the  pivot.  Notice  that  the  subarborescence  with 
root  c is  "rehung"  from  node  i,  and  that  nodes  on  the  pivot  stem  have  their 

predecessor  relationships  reversed.  The  flows,  X(  ) and  CPX(  ),  change  only  on 

the  backpaths  from  i and  j to  the  join.  The  simplex  multipliers,  U(  ),  and  the 

depth,  D(  ),  are  recomputed  for  all  nodes  in  the  subarborescence  with  root  c. 

The  most  expensive  part  of  the  pivot  is  the  update  of  the  simplex  multipliers  by 
the  addition  of  X to  the  U(  )'s  of  node  c and  all  nodes  in  the  subarborescence 
with  root  c.  With  the  data  structure  presented  so  far,  it  is  not  easy  to  identify 
all  the  (successor)  nodes  in  a subarborescence.  The  identification  of  these  nodes 
can  be  facilitated  by  a iA.avZKiCiZ  data  structure  that  begins  at  the  root  of  the 
predecessor  graph  and  exhaustively  "walks"  through  all  the  nodes  in  the  same  sequence 
that  the  nodes  occur  in  a triangulation  of  the  basis.  This  is  done  with  a node-length 
array  IT(  ) whose  i—  element  is  the  next  node  to  visit  from  node  i.  IT(  ) is 
thus  a different  way  to  represent  the  information  In  I.  It  is  convenient  to  make 

this  a circular  list  by  setting  the  IT(  ) of  the  last  node  in  the  triangulation 
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(Since  outgoing  arc  (3,4)  left  basis  at  its  capacity,  CP (2)  will  be  marked 
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Figure  4.  GNET  Arrays  and  Basis  after  Pivot  (Figure  2) 
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From  the  triangulation  construction  it  is  clear  that  a node  must  be  selected 
before  any  of  its  successors.  The  construction  of  the  triangulation  can  always  be 
modified  to  select  all  the  successors  (if  any)  of  a node  before  any  other  node  is 
considered.  For  this  restricted  class  of  triangulations  the  corresponding  traversal 
IT(  ) will  visit  all  the  nodes  of  a subarborescence  contiguously,  precisely  as  is 
required  in  the  update  of  the  simplex  multipliers.  This  traversal  is  the  obvious 
extension  to  the  m-ary  case  of  the  pfiZOUdzfL  (or,  equivalently,  dynastic  order) 
traversal  in  computer  science  literature.  The  recursive  definition  of  the  preorder 
traversal  given  in  the  computer  science  literature  [41]  reveals  precisely  its  value 
in  updating  the  simplex  multipliers. 

Let's  look  at  the  work  that  must  be  lone  in  the  pivot.  The  algorithm  visits 
each  node  of  the  subarborescence  with  root  c exactly  once.  It  proceeds  up  the 
pivot  stem  one  node  at  a time.  At  each  stem  node,  the  successors  of  the  next  lower 
stem  node  have  already  been  visited.  The  unvisited  successors  of  the  current  stem 
node  can  be  divided  into  two  groups:  the  nodes  visited  in  preorder  (by  iterating 

IT(  ))  from  the  stem  node  until  the  next  lower  stem  node  is  encountered,  called  the 
left  successors  of  the  stem  node,  and  the  remaining  unvisited  nodes  in  preorder  called 
the  right  successors  of  the  stem  node.  For  example,  in  Figure  2 the  left  successor  of 
stem  node  2 is  11,  and  the  right  successors  are  6 and  12.  Figure  5 gives  a general 
description  of  the  pivot. 

We  have  experimented  with  three  different  data  arrays  that  will  support  the 
onz-pam  pivot  computation:  depth  D(  ),  number  of  successors  NS(  ),  and  an 

additional  structure  PD(  ) to  be  discussed  below.  We  did  not  compare  another,  less 
efficient  two-pass  method  which  marks  nodes  with  sign  bits,  and  later  unmarks  them, 
much  like  the  procedure  described  for  locating  the  join  node.  These  arrays  are  used 
in  the  pivot  solely  to  answer  the  question  (Figure  5)  "An  unvisited  right  successor 
remains?";  this  question  is  nontrivial  because  the  node  information  now  available  to 
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US  Is  local  In  nature. 

Figure  6 shows  a detailed  description  of  the  pivot  traversal  using  depth.  The 
purpose  of  visiting  each  node  is  to  update  U(  ).  The  pivot  also  updates  the  arrays 
D(  ) and  IT(  ),  but  these  operations  are  omitted  for  clarity.  The  arrays  P(  ), 

X(  ) and  CPX(  ) are  recomputed  for  the  pivot  stem  only.  The  update  of  IT(  ) is 
easy  because  it  changes  only  for  the  nodes  on  the  pivot  stem  as  well  as  for  the  last 
left  successor  (if  any)  and  the  last  right  successor  (if  any)  of  each  pivot  stem  node. 
For  nodes  on  the  pivot  stem,  D(  ) is  updated  as  the  pivot  moves  up  the  stem.  The 
right  and  left  successors  of  each  pivot  stem  node  "inherit"  their  depth  from  the 
pivot  stem  node — if  the  depth  of  the  pivot  stem  node  changes  by  DADJ,  then  so  does 
the  dep_h  of  its  right  and  left  successors. 

The  number  of  successors,  NS(  ),  can  also  be  used  as  illustrated  in  Figure  7. 
The  updating  of  the  number  of  successors  is  not  shown  in  the  figure.  The  number  of 
successors  is  easy  to  update  because  it  changes  only  for  the  nodes  on  the  backpaths 
from  i and  j to  the  join.  Thus,  the  update  can  be  performed  for  stem  nodes 
simultaneously  within  the  pivot,  and  for  the  other  backpath  nodes  (from  i and  d) 
with  the  X(  ) and  CPX(  ) flow  update. 

Degeneracy  is  a critical  issue  in  transshipment  problems.  In  some  of  our  test 
problems  more  than  90%  of  the  (tens  of  thousands  of)  pivots  are  degenerate.  The 
search  for  the  join  may  be  aborted  when  degeneracy  is  encountered  since  the  only 
purpose  of  the  search  is  identification  of  the  leaving  arc  and  backpaths  for  flow 
change  (zero  in  the  case  of  degeneracy).  Stopping  short  also  tends  to  make  the  number 
of  nodes  visited  by  the  pivot  smaller.  In  the  depth  version,  the  lowest  degenerate 
arc  is  chosen  to  leave,  while  in  the  version  with  number  of  successors  the  smallest 
number  of  nodes  is  chosen  for  the  pivot.  But  the  successors  version  must  still 
locate  the  join  to  update  the  number  of  successors  on  each  backpath,  a relatively 
easy  process,  while  the  depth  version  requires  no  further  backpath  search. 
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Another  data  array  that  can  be  used  in  place  of  depth  or  number  of  successors 
is  position  in  the  triangulation,  or  equivalently  position  in  IT(  ).  For  each  node 
i define  'pKZOKdVl  dcitoncz  PD(i)  to  be  the  row  number  in  the  triangulation  of 

equation  i or  equivalently  to  be  the  number  of  iterations  of  IT(  ) beginning  with 

the  root  to  get  node  i.  Set  PD(root)  = 0. 

The  search  for  the  join  is  particularly  easy  with  preorder  distance.  The 
following  proposition  gives  a simple  construction  that  determines  the  join. 

Proposition  (McBride  [44],  Graves  and  Mcbride  [32]):  Given  a basis  with  a preorder 

triangulation,  for  any  two  nodes  i and  j such  that  PD(j)  > PD(i)  the  first  node 
h on  the  backpath  from  j with  PD(h)  ^ PD(i)  is  the  join. 

Proof:  The  construction  always  determines  a node  since  PD(root)  = 0.  The  join  is 

the  first  node  on  the  backpath  from  one  node,  say  j,  that  is  also  on  the  backpath 
from  the  other  node  i.  Thus  we  need  to  show  that  i is  either  h or  a sucessor 
of  h and  that  i is  not  a successor  of  any  other  node  on  the  backpath  from  j to 
h.  In  any  triangulation  PD(k)  > PD(P(k))  for  any  node  k,  thus  by  construction  h 
is  the  only  node  on  the  backpath  from  j to  h that  could  be  the  join.  In  a 

preorder  triangulation,  the  successors  of  any  node  k have  contiguous  PD(  ) numbers 

beginning  with  PD(k)  + 1.  Since  PD(h)  ^ PD(i)  < PD(j) , i either  equals  h or  is 
a successor  of  h.  Therefore,  h is  the  join  of  i and  j.  □ 

Using  PD(  ) the  outgoing  arc  can  be  determined  simultaneously  with  the  search 
for  the  join  and  if  degeneracy  is  discovered,  the  search  for  the  join  is  stopped. 

The  question  in  the  pivot  "An  unvisited  right  successor  remains?"  in  Figure  5 
is  equivalent  to  "PD(P(NT))  ^ PD(NR)?"  and  the  flowchart  for  preorder  distance  is 
equivalent  to  Figure  6 with  "D(NT)  > D(NR)"  replaced.  PD(  ) is  recomputed  for  all 
nodes  in  the  subarborescence  with  root  c and  some  additional  nodes. 

PRICING  MECHANISM 

The  pricing  operation  of  step  SI  in  the  simplex  algorithm  requires  a great  deal 
of  computational  effort.  As  in  other  large  scale  mathematical  programming  problems. 
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network  codes  can  spend  more  than  half  of  their  execution  time  selecting  incoming 
variables  by  pricing.  Thus,  the  pricing  mechanism  is  crucial  to  overall  performance 
and  must  also  provide  the  flexible,  broad  and  effective  external  controls  necessary 
to  permit  tuning  of  the  code  for  network  problems  with  special  or  even  bizarre  struc- 
ture. 

Choice  of  pricing  strategy  is  truly  an  art  for  large  scale  mathematical 
programming  codes.  It  is  based  on  intuition,  experience  and  empirical  experiments 
with  the  class  of  problems  to  be  solved.  Although  simplex  pricing  has  received  very 
little  exclusive  attention  in  the  literature,  fast  primal  codes  all  employ  some  form 
of  multiple  or  partial  pricing  of  subsets  of  the  variables  at  each  pivot  [47]. 

Examination  of  many  problem  trajectories  suggests  that  pricing  be  performed  in 
three  major  phases — an  opening  gambit,  a middle  game  and  an  end  game.  Initially,  a 
feasible  solution  is  being  constructed  and  pricing  must  select  incoming  variables 
carefully  among  many  available  choices  with  the  view  of  satisfying  violated  constraints. 

Later,  normal  gains  toward  optimality  seldom  justify  extensive  competition  among 
candidates.  Finally,  a favorable  candidate  becomes  quite  rare  and  thus  has  consider- 
able value  since  eventually  all  variables  mu.3t  be  exhaustively  priced  to  verify 
optimality. 

GNET  performs  pricing  by  selectively  using  a gznzfLdt  6can  mechanism  and  keeps  a 
set  of  good  pivot  information  in  a condidcutz  qtiZUZ.  For  a given  head  node,  the  scan 
selects  the  single  Incident  arc  (if  any)  pricing  most  favorably  (for  each  tail  examined 
incident  to  a particular  head  node  only  one  addition  and  one  comparison  are  required) 
and  places  it  on  the  candidate  queue.  The  candidate  queue  is  a varying  length  cyclic 
list  of  "interesting  node"  and  "good  arc"  information.  Each  entry  includes  a head 
node  and  either  a tail  pointer  specifying  the  location  of  a good  candidate  arc  or  a 
null  pointer  indicating  an  interesting  node  which  has  not  yet  been  priced.  The 
candidate  queue  mechanism  provides  for  user  control  of  network  solutions.  Although 
there  are  many  special  uses  for  these  controls,  for  example,  basis  crashing,  it  is 
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During  the  opening  gambit,  if  no  favorable  arc  is  discovered  after  examining  a 
block  on  the  candidate  queue,  another  block  of  entries  is  accessed,  and  so  forth, 
until  an  incoming  arc  is  discovered  or  the  candidate  queue  is  completely  emptied.  An 
exhausted  queue  is  refreshed  by  directing  the  general  scan  to  a specified  number 
("page"  IPG)  of  head  nodes. 

After  the  opening  gambit,  the  nodes  of  the  incoming  arc  are  no  longer  inserted 
in  the  candidate  queue.  After  each  complete  cycle  around  the  queue  it  is  refreshed 
by  a general  scan  of  IPG  nodes.  Each  of  these  scans  begins  where  the  previous  scan 
ended  and  thus  moves  cyclicly  through  the  arc  list  pricing  arcs  in  contiguous  storage 
locations.  When  managed  by  these  procedures,  the  candidate  queue  generally  grows 
during  and  just  after  the  opening  gambit  and  then  shrinks  finally  to  a few  good 
candidates.  The  end  game  is  played  by  concentrating  on  these  last  few  candidate  arcs. 

A major  portion  of  the  experimentation  with  sample  problems  has  been  devoted  to 
tuning  the  three  pricing  parameters  NNE,  NNS,  and  IPG  with  the  objective  of 
estimating  optimal  settings  as  simple  functions  of  solution  progress  and  easily  measured 
problem  characteristics  such  as  numbers  of  nodes,  arcs,  supply  and  demand  nodes, 
degree  of  capacitation , cost  range,  and  so  forth.  Also,  the  sensitivity  of  performance 
to  problem  structure  has  been  investigated. 

High  resolution  internal  computation  timing  is  often  difficult  to  carry  out  on 
contemporary  multiprogramming  computer  systems.  Some  published  result?’  are  highly 
suspect  due  to  unfortunate  oversights  in  timing  design.  The  computer  timing  routines 
often  take  longer  to  execute  than  the  network  algorithm  coding  segments  under  study. 
Therefore,  in  order  to  isolate  times  of  pricing  and  pivots,  a special  experiment  was 
designed  to  minimize  relative  timing  error.  Sample  problems  were  run  to  completion 
by  GNET  under  fixed  parameter  settings,  each  producing  a single  gross  solution  time 
and  number  of  pivots.  Next,  histories  of  the  entering  arcs  for  the  solutions  were  fed 
to  a specially  modified  code,  essentially  eliminating  the  pricing  mechanism;  these 
solution  times  represent  almost  exclusively  pivot  time.  Also,  nontiraing  runs  were 
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made  with  another  code  with  extensive  internal  data  collection  for  detailed  perform- 
ance analysis. 
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Pricing  schemes  were  tried  ranging  from  "first  negative"  (terrible)  to  "most 
negative"  (worse) . Optimal  intermediate  settings  of  NNE  and  IPG  lead  to  a nearly 
equal  distribution  of  time  between  pricing  and  pivoting.  As  the  length  of  the 
opening  gambit,  NNS,  is  increased,  solution  time  is  greatly  reduced.  However, 
beyond  some  point,  times  again  go  up  as  the  candidate  queue  becomes  clogged  with  bad 
nodes  inserted  by  the  pivots.  The  size  of  the  candidate  queue  swells,  and  then 
declines  as  predicted,  with  the  maximum  size  a complicated  function  of  all  three 
paraii;eters.  As  expected,  several  complete  cyclic  sweeps  of  the  entire  arc  list  in 
blocks  of  IPG  head  entries  always  occur  late  in  the  end  game,  but  the  candidate 
queue  focuses  attention  on  the  best  remaining  candidates.  The  total  number  of  pivots 
is  very  sensitive  to  NNE.  However,  total  solution  time  is  relatively  stable  as  NNE 
and  the  other  tuning  parameters  are  varied  from  reasonably  good  default  settings. 
Experiments  to  change  NNE  dynamically  during  solution  progress  have  not  improved 
performance.  In  fact,  many  dynamic  tuning  schemes  have  been  tried  without  much 
success . 

Hundreds  of  such  calibration  runs  have  been  performed  in  an  attempt  to 
determine  empirically  the  form  of  the  response  surface  of  execution  time  as  a func- 
tion of  the  pricing  parameters.  To  date,  the  best  general  automatic  default  settings 
, for  transshipment  problems  are 

f 
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; NNE  = 32  entries  examined  in  each  candidate  queue  block, 

NNS  = 3m/4  pivots  in  the  opening  gambit, 

IPG  = ra/10  head  nodes  priced  to  refresh  the  candidate  queue 
in  each  general  scan  block. 

: These  suggested  settings  are  surprisingly  robust  for  a wide  variety  of  problems.  How- 

ever, for  particular  classes  of  problems  sharing  uncommon  structure,  specific  tuning 
can  achieve  remarkable  further  improvements. 
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Degeneracy  Is  common  in  many  problems  (commonly  90  percent  of  the  pivots);  since 
degenerate  pivots  require  less  work  in  steps  S2  and  S3,  they  can  actually  be  a bonus 
rather  than  a worry.  The  pricing  mechanism  helps  obviate  the  need  for  any  additional 
methods  to  avoid  cycling  by  constantly  shuffling  the  cyclic  pricing  agenda.  (No  cycl- 
ing has  ever  been  encountered  in  our  experiments.)  Should  terminal  cycling  occur, 
the  exact  integer  solutions  will  certainly  lead  to  an  infinite  computational  loop. 

In  such  an  (unlikely)  case,  specification  of  slightly  modified  pricing  parameters 
(especially  NNE)  will  almost  certainly  avoid  the  cycle.  Formal  techniques  for  deal- 
ing efficiently  wit'.i  degeneracy  and  cycling  of  primal  simplex  network  codes  invite 
further  research. 

COMPUTATIONAL  EXPERIENCE 

The  family  of  GNET  codes  first  presented  at  the  Spring  1975  ORSA/TIMS  meeting 
represents  the  state  of  the  art  in  fast  large  scale  minimum  cost  network  flow  systems. 
Tables  1 and  6 report  solution  times  and  numbers  of  pivots  for  a set  of  standard  test 
problems  [using  NETGEN,  39]  that  have  also  been  solved  by  other  contemporary  codes. 

These  solution  times  are  achieved  with  the  pricing  parameters  set  at  their  default 
values  as  in  the  code  that  is  distributed.  Thus,  our  experiments  are  completely 
reproducible  by  other  researchers.  Benchmarks  on  various  machines  (see  Table  2,  for 
instance)  generally  agree  with  standard  hardware  comparisons  of  computer  performance 
and  show  that  the  times  in  Tables  1 and  6 are  superior  to  all  published  times  of  which 
we  are  aware.  However,  we  cannot  verify  published  claims  of  machine,  compiler  and 
installation  performance  variations.  In  any  case,  future  codes  produced  by  incorporating 
fresh  research  ideas  will  undoubtedly  make  our  current  performance  records  obsolete. 

We  do  not  hold  that  primal  network  codes  are  anywhere  near  an  asymptotic  efficiency 
level . 

The  GNET  family  presently  includes  GNET/Depth,  GNET /Successors  and  GNET/Preorder- 
distance  as  described  above,  as  well  as  other  variants  to  be  discussed  below.  The 
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TABLE  1 
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GNET/Depth  Performance  on  Several  Transshipment  Examples 

(Default  Tuning) 

Percent  IBM  360/67 


Problem 

Nodes 

Sources 

Sinks 

Arcs 

Gap'd 

Seconds 

Pivots 

NG27 

400 

4 

12 

2,676 

80 

2.6 

607 

NG36 

8,000 

200 

1,000 

15,000 

0 

212 

13,012 

NG37 

5,000 

150 

800 

23,000 

0 

138 

11,610 

NG38 

3,000 

125 

500 

35,000 

0 

97 

10,637 

NG39 

5,000 

180 

700 

15,000 

0.7 

113 

9,553 

NG40 

3,000 

100 

300 

23,000 

0.7 

67 

6,409 

NPS 

10,000 

50 

5,000 

21,000 

100 

441 

22,153 

XNl 

5,000 

100 

4,800 

40,000 

100 

290 

12,111 

TABLE  2 

GNET/Depth  Performance  on  'NG27' 

- Calibration  with  a 

Highly  Capacitated  Transshipment  Problem 

Machine 

Solution  Seconds 

CDC  7600 

0.3 

IBM  370/168 

0.5 

CDC  6600 

0.6 

Ti  ASC 

0.7 

IBM  360/91 

0.8 

UNI  1108 

2.2 

IBM  360/67 

2.3 
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three  basic  versions  of  GNET  have  been  tested  on  a suite  of  randomly  generated  problems 
and  real  formulations  to  determine  which  version  is  best.  Experiments  on  smaller  pro- 
blems (less  than  500  nodes  and  5,000  arcs)  show  that  all  three  are  remarkably  close 
with  the  depth  version  a narrow  favorite  for  sheer  speed.  Although  preorder  distance 
is  preferable  for  mathematical  reasons  when  extending  the  data  structures  to  non- 
network problems,  large  scale  testing  has  been  limited  to  depth  and  number  of  succesors 
versions.  GNET/Successors  execution  times  differ  from  the  GNET/Depth  times  by  at  most 
2 percent  on  the  problems  in  Table  1. 

All  versions  begin  with  a read  and  edit  routine  and  an  arc  list  sort  to  create 
the  compressed  arc  arrays.  Supplies  and  demands  are  determined  for  each  node  and  the 
candidate  queue  is  loaded  with  the  demand  nodes.  Simplex  pricing  via  the  candidate 
queue  mechanism  is  used  to  reduce  artificial  flows  to  zero  rather  than  a Phase  I- 
Phase  II  approach  (used  in  earlier  GNET  versions).  The  high  cost  for  artificial  arcs 
is  computed  as  the  product  of  the  number  of  nodes  (maximum  path  length)  and  the 
maximum  absolute  arc  cost,  thus  guaranteeing  that  a feasible  problem  will  have  no  flow 
on  artificial  arcs  in  a final  basis.  The  cost  is  attached  only  to  artificial  arcs 
incident  to  demand  nodes;  experiments  attaching  the  cost  to  arcs  incident  to  supply 
nodes  alone,  and  to  both  supplies  and  demands  have  not  improved  performance.  GNET  is 
tuned  for  node  numbering  ascending  from  supply  through  transshipment  to  demand  nodes, 
and  departures  may  somewhat  degrade  performance  of  the  demand  driven  pricing  mechanism 
(for  general  transshipment  problems  the  general  scan  blocks  are  contiguous  and  exhaus- 
tive). Another  underlying  assumption  has  been  that  problems  will  normally  have  many 
more  sinks  than  sources.  Problems  with  more  sources  than  sinks  or  erratic  node 
numbering  can  usually  be  easily  reformulated  to  our  preference  if  solution  speed  is 
of  prime  importance;  with  large  scale  models  this  is  a minor  undertaking  in  the  pro- 
blem generator  software. 

Use  of  random  test  problems  in  tuning  network  codes  has  its  pitfalls.  GNET 
solves  real  network  models  much  faster  than  random  NETGEN  problems  of  nominally 


comparable  size  and  structure;  this  suggests  that  much  remains  to  be  learned  from 
further  investigation  of  special  problem  structure.  Large  random  test  problems  are 
very  expensive  to  generate  with  NETGEN,  requiring  about  five  times  more  computer  time 
than  the  associated  C2^ET  optimization  (but  much  less  region).  However,  the  reproduci- 
bility of  the  experiments  and  results  is  so  important  that  it  justifies  the  expense. 
Also,  the  cost  and  awkwardness  of  using  magnetic  tapes  with  standard  problem 
libraries  is  avoided.  Unfortunately,  large  NETGEN  problems  can  vary  slightly  between 
computers  with  the  length  of  the  real  mantissa  (this  is  due  to  a random  number  gene- 
rator which  simulates  35  bit  integer  arithmetic,  normalization  of  the  integer  result 
to  a real  0-1  variable  on  the  host  machine  and  subsequent  transformation  back  to  a 
uniform  integer  with  the  desired  range).  NETGEN  puts  capacities  of  +1  on  all  arcs  in 
assignment  problems,  needlessly  complicating  their  solution  by  a general  capacitated 
network  code.  Negative  demands  can  also  be  generated.  These  and  other  minor  NETGEN 
problems  arc  overcome  by  a few  program  modifications.  It  also  would  be  worthwhile  to 
add  the  capability  to  generate  multiple  arcs  and  special  problem  structures. 

GNET  execution  times  do  not  show  much  sensitivity  to  cost  ranges  contrary  to 
past  reports  [55].  Experiments  were  performed  in  which  individual  problems  were 
solved  repeatedly  with  only  the  costs  modified  either  by  low  order  digit  truncation 
or  addition  of  uniform  random  low  order  digits  as  necessary  to  provide  the  desired 
cost  range.  Solution  times  are  very  insensitive  to  cost  ranges  thus  produced,  seldom 
varying  by  more  than  15  percent.  The  outcome  of  one  such  expt^riment  is  shown  in 
Table  3.  This  surprising  result  appears  to  be  attributable  to  the  candidate  queue 
pricing  mechanism. 

Uncapacitated  transshipment  problems  are  input  to  GNET  with  upper  bounds  on  each 
arc  exceeding  total  problem  supplies.  Minor  modifications  of  GNET  can  reduce  solution 
times  by  about  10  percent  for  strictly  uncapacitated  problems;  a more  important  issue 
is  the  potential  elimination  of  the  arc-length  array  of  capacities.  This  space 
economy  may  also  be  realized  on  lightly  capacitated  problems  by  modification  of  the 
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GNET/Depth  Performance 

TABLE  3 

for  'NG40'  with  Varied  Ranges 

Random  Uniform  Costs  (Original 

Range  1-100) 

i 

Cost  Range 

IBM  360/67 
Seconds 

Pivots 

1-10 

84 

6,726* 

1-100 

76 

6,933 

c 

1-1,000 

74 

6,921^ 

1-10,000 

74 

6,958^ 

Low  order  digit  truncated  from  original  problem  costs. 


Low  order  digits  generated  by  RANDU  and  concatenated  with 
original  problem  costs. 


TABLE  4 


Lightly  Capacitated  Transshipment  Problems  - 
Performance  of  GNET/Depth  and  Modifications 
Solving  the  Equivalent,  Reformulated  uncapacitated  Problems  [11] 


NT39 

o 

•-> 

. 

IBM  360/67 
Seconds 

Pivots 

IBM  360/67 
Seconds 

Pivots 

GNET (depth) 

113 

9,553 

67 

6,409 

Preprocess 

113 

9,615 

71 

7,511 

Transform  on-the-fly 

114 

9,652 

65 

6,733 

'f' 
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capacity  array  or  by  problem  reformulation.  One  section  of  a report  by  Cheong  [11] 
includes  a catalog  of  problem  transformations  useful  in  network  models  and  gives 
results  of  experiments  with  several  large,  lightly  capacitated  transshipment  problems. 
Two  special  versions  of  GNET/Depth  were  prepared  which  utilize  a well-known  transfor- 
mation to  replace  each  capacitated  arc  by  a pair  of  uncapacitated  arcs  and  a new  node. 
(The  tail  of  the  replaced  arc  is  moved  to  the  new  node  and  supply  equal  to  the  replaced 
capacity  is  also  shifted.  The  new  node  is  then  connected  to  the  old  tail  node  by  an 
artificial  arc  with  zero  cost.)  One  modification  performs  the  transformation  to  the 
arc  list  before  solution.  The  other  carries  out  the  transformation  "on-the-fly"  as 
arcs  with  capacities  are  introduced  into  the  basis.  Table  4 gives  an  example  of  per- 
formance for  the  codes,  with  the  modifications  each  using  one  third  less  space  for 
the  arc  arrays  with  little  speed  degradation.  As  the  proportion  of  capacitated  arcs 
increases,  the  space-time  tradeoff  becomes  much  less  favorable. 

POSTOPTIMAL  ANALYSIS,  REOPTIMIZATION  AND  FURTHER  REFINEMENTS 

In  some  applications,  analysis  of  the  sensitivity  of  the  optimal  solution  to 
modifications  in  the  supplies,  demands  and  cost  coefficients  is  desirable.  All  the 
postoptimal  analysis  for  linear  programming  problems  can  be  done  with  primal  trans- 
shipment codes.  The  data  structures  described  above  support  efficient  techniques  for 
this  analysis.  Ranging  [47]  of  problem  coefficients  traditionally  requires  informa- 
tion from  the  inverse  of  the  optimal  basis;  columns  are  required  for  ranging  of 
supplies  and  demands,  rows  are  required  for  ranging  of  costs  of  basic  arcs.  As 
described  above,  the  predecessor  function  can  generate  columns  of  the  inverse  and 
the  traversal  mechanism  can  generate  rows.  Note  that  while  the  ranging  of  cost 
coefficients  of  basic  arcs  is  simple  in  principle  and  much  faster  than  for  linear 
programming  systems,  for  applications  with  many  arcs  it  can  be  time  consuming  rela- 
tive to  the  time  to  construct  the  optimal  solution. 
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A more  important  issue,  especially  for  applications  that  have  a network  embedded 
in  a larger  problem,  is  the  reoptimization  of  a problem  after  modifications  to  the 
problem  coefficients.  The  case  for  primal-dual  or  dual  algorithms  is  often  based  on 
the  reputed  ease  of  reoptimization;  although  reoptimization  is  conceptually  easier 
than  for  primal  algorithms , this  does  not  imply  that  primal-dual  or  dual  algorithms 
reoptimize  more  quickly.  There  have  not  been  comparisons  of  primal  and  nonprimal 
algorithms  for  reoptimization  of  large  scale  problems  (indeed,  it  is  not  clear  what 
types  of  problems  would  constitute  a fair  comparison).  In  any  case,  the  computer 
storage  advantage  of  primal  algorithms  is  still  maintained  since  reoptimization 
requires  no  additional  data  arrays. 

The  following  design  considerations  were  used  for  a GNET  reoptimization  procedure: 
1)  there  would  be  no  modifications  to  the  primal  optimization  segments  of  GNET,  2)  no 
new  arcs  would  be  inserted  into  the  arc  list,  and  3)  the  initial  basis  for  the  reopti- 
mization would  include  as  much  of  the  previous  optimal  basis  as  possible.  Algebraic- 
ally, the  changes  to  the  supplies,  demands,  lower  bounds  and  upper  bounds  are 
translated  into  a vector  of  changes,  d,  to  be  added  to  the  original  right  hand  side 
vector  b.  Using  the  previous  optimal  basis  B,  the  set  of  equations  Bk  = b + d is 
solved  with  artificial  arcs  introduced  into  the  basis  as  necessary  to  avoid  an 
infeasible  x.  Since  the  previous  optimal  solution  x is  available  (b  is  not) , 

Bx  = d is  solved  and  then  added  to  x to  get  k. 

The  vector  d originates  from  coefficient  modifications  by:  1)  changes  in 
supplies  and  demands,  2)  a change  of  A in  the  upper  bound  of  an  arc  (i,j)  out  of 


the  basis  at  this  upper  bound  (A  is  added  to  d^  and  subtracted  from  ^ 


change  of  A in  the  lower  bound  of  an  arc  (i,j)  out  of  the  basis  at  zero  (A  is 


subtracted  from  d^  and  added  to  ^ traversal  mechanism  for  back  substitution 


of  the  basis  is  provided  by  a single  pass  through  the  IT(  ) array.  The  vector  x 
is  constructed  by  backsolving  with  the  predecessor  array  P(  ) , at  each  step  k^  is 


computed  as  x^  + x^.  If  k^  is  nonnegative  and  less  than  or  equal  to  its  capacity. 
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the  backsolvlng  continues.  Otherwise,  an  artificial  vector  from  the  node  to  the  root 
Is  Introduced  to  replace  the  arc  which  leaves  the  basis  at  zero  (If  < 0)  or  Its 
capacity  (If  > capacity) . The  artificial  arc  Is  oriented  so  that  Its  flow  Is 
positive  and  It  Is  assigned  the  large  cost.  The  exchange  of  an  arc  and  an  artificial 
arc  Involves  "rehanglng"  the  node  and  all  Its  successors;  this  Is  accomplished  with 
only  three  changes  In  the  preorder  traversal  IT(  ) by  Inserting  these  nodes  at  the 
end  of  IT(  ).  Arcs  removed  from  the  basis  are  placed  In  the  candidate  queue.  The 
backsolvlng  continues  until  the  whole  basis  Is  recomputed.  If  any  artificial  arcs 
have  been  Introduced  it  is  necessary  to  recompute  the  simplex  multipliers.  The  pro- 
blem is  then  reoptimized  by  GNET. 

If  cost  coefficient  changes  Involve  only  nonbaslc  arcs,  reoptimization  Is 
necessary  only  if  one  or  more  now  price  favorably  (such  arcs  are  added  to  the  candi- 
date queue).  If  cost  coefficients  change  for  basic  arcs  it  is  necessary  to  recompute 
the  simplex  multipliers  and  then  reoptimize  with  GNET. 

Since  GNET  is  so  fast,  experience  Indicates  that  if  there  are  many  coefficient 
changes  it  may  be  more  efficient  to  begin  with  an  all  artificial  basis  with  the 
previous  optimal  basis  arcs  preloaded  In  the  candidate  queue.  On  the  other  hand, 
applications  often  require  many  particular,  recurring  reoptimizations  of  some  special 
type  which  permit  more  efficient  (and  problem  specific)  methods  to  be  applied. 

One  of  the  major  Ideas  of  mathematical  programming  Is  that  elements  of  the 
simplex  tableau  may  be  generated  as  needed  rather  than  stored  and  updated  at  each 
pivot.  For  example,  as  described  above,  in  GNET  the  columns  of  the  basis  inverse  are 
generated  when  needed  to  determine  the  outgoing  arc  by  iterating  the  predecessor 
function.  However,  the  basic  flows  and  simplex  multipliers  are  explicitly  stored  and 


the  ones  that  change  are  updated  during  each  pivot.  Analysis  of  computational  experi- 
ments show  that  the  major  portion  of  the  calculations  in  the  pivot  (step  S3)  is  to 
update  the  simplex  multipliers  U(  ),  depth  D(  ) and  the  preorder  traversal  array 


IT  ( ) (note  that  IT(  ) is  maintained  solely  to  allow  efficient  update  of  the 
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multipliers) . The  authors  were  thus  led  to  consider  modifications  to  GNET  that  would 
allow  some  or  all  of  the  multipliers  to  be  generated  as  needed. 

It  follows  immediately  from  (11)  that  for  basis  arc  k from  node  i to  j: 

C(k)  - U(i)  + U(j)  = 0.  Arbitrarily  setting  U(root)  = 0,  it  is  possible  to  solve 
for  all  the  U(  )'s  by  front  substitution  with  the  triangulated  basis;  that  is, 
solve  for  the  U(  ) in  the  order  IT(root) ,IT(IT(root)) , etc.  To  compute  a particular 
U(h)  it  is  necessary  only  to  calculate  some  of  the  other  U(  )'s,  namely  those  for 
the  nodes  on  the  backpath  from  h to  the  root.  If  k is  the  arc  joining  node  h to 
its  predecessor  P(h),  then  U(h)  = U(P(h))  ± C(k)  with  "+"  if  the  arc  is  oriented 
from  h to  P(h)  and  otherwise.  Since  there  are  many  nodes  in  large  problems 

and  only  a relatively  few  multipliers  change  from  pivot  to  pivot,  it  does  not  seem 
worthwhile  to  generate  all  multipliers  at  each  pivot.  Rather,  we  choose  to  store 
explicitly  enough  of  the  multipliers  so  that  for  each  node  i either  U(i)  is  stored 
explicitly  or  U(P(i))  is  stored  explicitly.  In  the  latter  case  a single  addition 
generates  the  multiplier  only  when  it  is  needed. 

The  impetus  for  this  approach  comes  from  consideration  of  the  capacitated  trans- 
portation problem  with  many  demand  nodes  (sinks)  and  relatively  few  supply  nodes 
(sources).  This  is  an  important  special  case  of  the  minimum  cost  network  model  that 
has  many  applications.  A typical  model  is  a distribution  system  with  a few  plants  and 
many  warehouses  or  a few  centralized  warehouses  and  many  customers.  This  is  the  type 
of  problem  most  often  encountered  in  multicommodity  distribution  problems,  e.g.  [22]. 
Another  application  that  requires  the  repeated  solution  of  many  such  problems  is  the 
traffic  assignment  problem. 

A special  version  of  GNET/Depth  called  TNET  was  developed  for  this  problem.  The 
multipliers  are  explicitly  kept  only  for  the  relatively  few  sources  of  the  problem. 
Since  in  the  (bipartite)  transportation  problem  arcs  only  join  sources  to  sinks,  the 
predecessor  of  each  sink  node  must  be  a source  and  thus  the  multiplier  for  each  sink 
is  computed  by  a single  addition  (also,  the  pricing  mechanism  loads  only  sinks  as 
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interesting  nodes  during  the  opening  gambit).  The  GNET  arc  storage  is  ideal  for  this 
calculation  since  all  the  arcs  with  sink  node  j are  stored  contiguously;  thus  in  the 
general  scan  all  are  priced  out  together;  and  only  a single  calculation  of  the  sink 
multiplier  is  required. 

Since  the  traversal  array  1T{  ) is  maintained  solely  for  the  purpose  of 
updating  the  multipliers,  IT(  ) and  D(  ) need  to  be  maintained  only  for  the 
relatively  few  source  nodes.  It  is  also  convenient  to  include  in  IT(  ) the  rela- 

1 

tively  few  sinks  that  have  a successor.  An  alternate  description  is  that  IT(  ) is 
maintained  only  for  the  smallest  possible  subarborescence  that  includes  all  the 
sources;  this  is  easy  to  maintain  since  at  each  pivot  the  only  possible  nodes  that 

^ can  join  the  subarborescence  are  the  two  ends  of  the  incoming  arc  and  the  only 

[ possible  nodes  to  leave  are  the  two  ends  of  the  outgoing  arc.  No  additional  storage 

[ is  required,  the  U(i)  for  sink  i is  used  to  store  the  cost  of  the  arc  from  P(i) 

I. 

j 

j to  i and  the  sinks  not  in  the  subarborescence  are  marked  with  0 in  D(  ) . The 

i 

I predecessor  array  P(  ) is  maintained  exactly  as  before  and  the  determination  of  the 

i 

t backpaths  and  outgoing  arc  is  the  same. 

Experiments  (see  Table  5)  show  that  TNET  is  significantly  faster  than 
GNET  for  transportation  problems  with  many  more  sinks  than  sources.  The  pivot 
choice  is  the  same  in  both  the  original  and  modified  versions  so  with  the  same 

pivot  choice  is  the  same  in  both  the  original  and  modified  versions  so  with  the  same 
value  of  the  pricing  parameters  the  same  sequence  of  pivots  is  generated.  For 
NNE  = 32,  the  reduction  in  time  for  TNET  is  due  solely  to  savings  in  the  update  of 
U(  ),  IT(  ),  and  D(  ) in  step  S3  (step  S2  is  identical  and  step  SI  is  slightly 
slower  for  TNET).  Experiments  with  six  smaller  problems  show  that  GNET  spends  roughly 
half  its  time  in  step  SI;  with  TNET  there  was  80%  less  time  for  S3  and  5%  more  time  in 
SI.  Since  the  pivot  step  S3  is  so  much  faster  in  the  modification  it  is  not  worth- 
while to  select  the  incoming  arc  as  carefully  as  in  GNET;  experiments  show  that 
NNE  =8  is  a good  setting  for  TNET.  Table  5 shows  that  with  NNE  = 8 there  are 
significantly  more  pivots  than  with  NNE  = 32,  but  total  solution  time  is  reduced. 
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TABLE  5 

Transportation  Problems  with  Relatively  Few  Sources 


Problem 

Sources 

Sinks 

Arcs 

NNE 

IBM  360/67  Seconds 
TNET  GNET 

Pivots 

TN7 

100 

4,900 

40,000 

32 

127 

274 

11,909 

8 

99 

17,583 

TN8 

250 

4,750 

40,000 

32 

142 

285 

12,910 

8 

135 

19,759 

TN9 

10 

4,990 

40,000 

32 

104 

261 

9,574 

8 

79 

12,327 

Generated 

with  NETGEN  [39], 

all  arcs 

are 

capacitated, 

cost  ranges  are  0-100  (TN7)  and  0-1000  (TN8,TN9). 


TABLE  6 

XNET/Depth  Performance  with  NNE  = 8 


Problem 

IBM  360/67 
Seconds 

Pivots 

Problem 

IBM  360/67 
Seconds 

Pivots 

NC27 

2,7 

964 

NPS 

265 

29,045 

NG36 

111 

17,993 

XNl 

136 

19,726 

NG37 

110 

18,195 

NG38 

94 

13,124 

TN7 

112 

17,583 

NG39 

76 

14,809 

TN8 

148 

19,759 

NG40 

59 

11,002 

TN9 

92 

12,327 
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Although  TNET  is  designed  for  transportation  problems  with  many  sinks  and 
relatively  few  sources,  it  does  well  (about  the  same  as  GNET)  on  problems  with  an 
equal  number  of  sources  and  sinks.  TNET  handles  problems  with  many  sources  and  few 
sinks  by  a minor  modification  in  the  input  that  reverses  the  orientation  of  the  arcs. 

The  same  idea  is  extended  to  the  general  minimum  cost  network  problem  in  a pro- 
gram called  XNET.  The  multipliers  are  explicitly  kept  only  for  nodes  that  have 
successors  in  the  predecessor  graph.  For  example,  in  Figure  2 IT(  ) is  maintained 
only  for  nodes  13,  8,  1,  4,  3,  2,  6 in  that  order.  The  U(  ) and  D(  ) arrays  for 
nodes  with  no  successors  are  used  as  in  TNET.  In  order  to  keep  the  multipliers  only 
for  the  nodes  with  successors,  it  is  necessary  to  maintain  an  array  that  records  for 
a node  the  number  of  its  successors  that  do  not  have  explicit  multipliers  ("aggregated 
successors"  A(  )). 

The  results  in  Table  6 indicate  that  XNET  is  faster  than  GNET  on  all  problems 
and  significantly  faster  on  problems  with  relatively  many  sinks.  XNET  is  only 
slightly  slower  than  TNET  on  the  problems  that  TNET  is  tailored  for,  but  the  loss  on 
these  problems  is  balanced  by  its  generality  and  by  its  dominance  of  GNET.  XNET  is 
successful  because  the  predecessor  graphs  have  many  nodes  with  no  successors.  In  many 
practical  problems  known  to  the  authors,  most  nodes  are  pure  sinks;  since  the 
successor  in  the  predecessor  graph  of  a pure  sink  must  be  a node  that  is  not  a pure 
sink,  many  pure  sink  nodes  have  no  successors. 

TNET  and  XNET  are  refinements  of  GNET/Depth;  number  of  successors  or  preorder 
distance  could  also  be  used  for  such  modifications. 

A modified  version  of  XNET  is  also  being  used  on  a set  of  large,  100-percent 
dense  problems.  The  network  models  are  uncapacitated  single  commodity  transportation 
problems  embedded  in  a recent  implementation  of  a multicommodity,  multiple  time  period 
econometric  model  described  in  [34].  A prototype  problem  size  is  200  sources  and  300 
sinks  t and  thus  60,000  arcs)  with  a Mtrix  of  region  to  region  bulk  transport  costs. 
The  XNET  modificaticn  is  stripped  of  arc-length  arrays,  list  references  are  modified 
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for  the  cost  matrix  and  several  capacitated  features  removed.  Also,  the  candidate 
queue  is  modified  to  access  only  sink  nodes  and  to  price  a restricted  menu  of  the 
few  (ultimately  5)  cheapest  sources  for  each  sink  during  the  opening  gambit.  Optimi- 
zation from  a cold  start  requires  8.6  seconds  (2350  pivots)  on  an  IBM  370/168. 

Tuning  of  pricing  parameters  produces  surprisingly  little  improvement.  Reoptimiza- 
tion procedures  are  employed  to  exploit  period-to-period  similarity  of  optimal  bases 
(despite  significant  temporal  variations  in  problem  structure).  A typical  reoptimi- 
zation time  from  a crashed  basis  is  1.2  seconds. 

HISTORICAL  PERSPECTIVE 

Although  the  authors  have  not  had  access  to  any  other  primal  transshipment 
computer  programs,  it  is  possible  to  identify  some  of  the  major  ideas  from  papers  and 
presentations  [10,  26,  27,  28,  29,  32,  45,  46,  54,  55].  The  major  design  decisions 
and  important  coding  specifics  may  vary  widely  in  these  systems,  but  all  the  success- 
ful contemporary  large  scale  codes  seem  to  be  based  on  a few  key  ideas.  Notwith- 
standing our  limited  knowledge  of  specific  techniques  used  by  others,  we  trace  the 
development  of  these  ideas  as  best  we  can. 

The  major  ideas  underlying  all  contemporary  primal  network  codes  are  the 
representation  of  the  basis  as  a predecessor  graph,  a traversal  mechanism  to  update 
simplex  multipliers,  the  use  of  depth  or  number  of  successors  or  preorder  distance  to 
determine  the  outgoing  arc  and  to  facilitate  the  update  of  the  simplex  multipliers, 
the  use  of  a pricing  mechanism  such  as  the  candidate  queue,  and  generation  of 
simplex  multipliers. 

The  predecessor  array  is  used  in  all  the  codes.  This  construction  in  primal 
network  codes  goes  back  at  least  to  Glicksman,  L.  Johnson  and  Eselson  [25]  1960  and 
goes  back  further  in  other  disciplines.  It  is  a standard  approach  in  many  research 
areas . 
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Srinivasan  and  Thompson  [54J  1972  have  proposed  the  use  of  depth  to  identify  the 
backpaths  and  to  determine  the  outgoing  arc.  This  has  been  adopted  by  Glover,  Karney 
and  Klingman  126]  1974,  Mulvey  145]  1974,  [46]  1975  and  GNET/Depth  [6]  1975.  The 
method  [54]  for  identifying  the  join  of  the  backpaths  is  exactly  as  described  for 
GNET,  however,  the  update  of  depth  at  each  pivot  in  that  paper  is  a more  involved 
approach  that  does  not  aid  in  the  pivot  as  described  here. 

Preorder  traversal  was  first  used  in  transshipment  codes  by  Glover,  Klingman  and 
Stutz  [29]  1974  who  called  it  the  Augmented  Threaded  Index  method  (ATI).  They  show 
that  it  is  more  efficient  than  the  triple  labeling  scheme  of  E.  Johnson  [37].  Indepen- 
dently, McBride  [44]  1973  and  Graves  and  McBride  [32]  1973  have  developed  the  zero  to 
the  right  (ZTR)  triangulation  of  network  bases  which  has  since  been  shown  to  be 
equivalent  to  preorder.  Preorder  has  been  used  in  the  computer  science  literature 
for  at  least  20  years  and  has  been  used  in  contemporary  shortest  path  algorithms  [24] 
1973.  Preorder  has  been  adopted  by  Srinivasan  and  Thompson  and  Mulvey  and  it  is  used 
in  all  versions  of  GNET.  As  discussed  above  and  shown  in  Figure  6,  depth  makes 
possible  a one-pass  update  of  the  simplex  multipliers  that  simultaneously  updates  the 
preorder  and  depth  (or  number  of  successors  or  preorder  distance) . Since  the  pivot 
taxes  much  more  time  than  the  determination  of  the  outgoing  arc,  the  use  of  depth  in 
the  update  of  the  simplex  multipliers  is  more  critical  to  the  success  of  GNET  than  its 
use  to  find  the  outgoing  arc.  We  do  not  know  the  pivot  details  or  the  use  of  depth 
(if  any)  in  other  codes  nor  do  we  know  if  any  have  a one-pass  update. 

Srinivasan  and  Thompson  [54]  1972  have  proposed  number  of  successors  to  be  used 
with  depth  to  determine  how  many  nodes  are  in  the  subarborescence  that  is  rehung  in 
the  pivot.  If  the  subarborescence  has  less  than  half  the  nodes  of  the  problem,  they 
propose  moving  the  root  and  performing  the  pivot  or.  the  smaller  part  of  the  prede- 
cessor graph.  Our  experiments  show  that  the  subarborescence  can  always  be  expected  to 
have  significantly  less  than  half  the  nodes,  so  we  have  rejected  the  idea  of  moving 
the  root.  However,  we  have  discovered  that  number  of  suu  :essors  can  be  used  in  a way 


not  evisaged  by  Srinivasan  and  Thompson.  As  described  above  it  can  replace  depth  to 
support  the  determination  of  the  outgoing  arc  and  the  one-pass  update  of  the  simplex 
multipliers.  We  have  also  noted  that  the  number  of  successors  for  nodes  in  the  pivot 
subarborescence  changes  only  for  nodes  on  the  pivot  stem  avoiding  the  update  indicated 
in  [55].  Subsequent  to  our  presentation  of  GNET/Successors , Glover  and  Klingman  [28] 
1975  have  proposed  a number -of -successors  version  of  their  algorithm.  They  describe 
a one-pass  update  of  the  simplex  multipliers  that  uses  an  additional  node  length 
array  to  store  for  each  node  the  last  successor  ranked  by  preorder.  They  conjec- 
ture [28,  p.  7]  that  a code  based  on  this  approach  will  dominate  their  previous  codes. 
This  may  be  true  for  their  codes,  but  as  discussed  above  this  is  not  our  experience 
in  comparing  our  one-pass  depth  version  against  our  one-pass  mamber  of  successors 
version.  Also, as  indicated  in  Figure  7, GNET/Successors  accomplishes  a one-pass  update 
without  an  additional  array  to  store  and  manipulate. 

Preorder  distance  was  used  by  McBride  [44]  1973  and  Graves  and  McBride  [32]  1973. 
They  developed  the  zero-to-right  property  and  established  the  proposition  stated  above 
to  find  the  join. 

Mulvey  [45]  uses  a candidate  list  of  arcs  that  corresponds  to  the  "partial 
suboptimization"  that  is  used  in  commercial  linear  programming  systems  [47].  His 
I candidate  list  is  controlled  by  two  parameters:  and  x^.  In  our  terminology, 

: general  scans  of  nodes  are  done  until  there  are  arcs  on  the  candidate  list, 

then  up  to  x^  pivots  are  performed  by  choosing  arcs  from  the  candidate  list.  The 
candidate  list  is  then  discarded  and  the  process  begins  again.  As  noted  in  [10], 
Mulvey 's  approach  has  been  adopted  by  Glover  and  Klingman. 

The  candidate  queue  used  in  GNET  has  evolved  from  similar  mechanisms  developed 
by  Graves  for  more  general  mathematical  programming  systems.  It  is  unique  in  that  it 

I contains  both  nodes  and  arcs,  it  is  used  cyclically  and  arcs  are  never  removed  as  long 

as  they  continue  to  price  out  favorably. 
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The  Srinivasan  and  Thompson  code  [55J  and  the  Harris  code  [33J  are  for  dense 

uncapacitated  transportation  problems;  all  other  contemporary  codes  known  to  the  authors 
solve  the  capacitated  transshipment  problem  with  a sparse  representation  of  the  net- 
work similar  to  that  described  here. 

Harris  [33]  1976  has  described  a code  for  dense  uncapacitated  transportation 
problems  that  have  few  sources  and  many  sinks.  The  brief  description  indicates  that 
the  simplex  multipliers  are  not  stored  for  the  sinks  but  there  are  no  details  on  the 
handling  of  the  pivot  nor  on  the  use  of  a traversal  mechanism.  Independently,  the 
authors  have  developed  TNET  for  the  sparse  capacitated  transportation  problem  and 
subsequently  developed  XNET  for  the  general  capacitated  transshipment  problem.  TNET 
and  XNET  are  specific  examples  of  the  generation  of  simplex  multipliers;  other 
refinements  that  generate  more  (or  all)  the  multipliers  are  possible. 

There  is  considerable  literature  on  postoptimal  analysis  and  parametric 
programming  for  transportation  and  transshipment  problems  (e.g.  [9,  52,  53]). 

CONCLUSION 

The  GNET  programs  are  small,  fast  and  easy  to  modify.  They  integrate  many  well- 
known  techniques  from  mathematical  optimization  and  computer  science.  It  is  important, 
however,  to  carefully  discriminate  between  the  underlying  ideas  of  mathematical  program- 
ming and  the  computer  science  topics  applied  to  their  implementation.  A sound 
mathematical  footing  is  required  for  generalization  beyond  pure  networks.  The 
computer  science  literature  has  contributed  a great  deal  to  the  local  representation 
of  global  information  in  trees  and  graphs,  and  has  given  valuable  recursive  methods  for 
manipulating  data  structures;  many  of  these  techniques  can  be  applied  to  the  basic 
arborescence.  But  some,  such  as  rerooting,  pruning,  balancing  and  even  conversion  to 
an  equivalent  binary  tree  with  extra  dummy  nodes  and  arcs,  do  not  work  well  for  network 
problems.  In  these  cases,  the  mathematical  interpretation  implies  (as  does  the  experi- 
mental evidence)  that  these  devices  are  needless  complications  that  increase  solution 
expense  or  introduce  superfluous  equations  and  variables. 
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In  our  experience  large  scale  problems  are  always  created  from  source  data  by 
problem  generator  programs.  The  problem  generator  may  occupy  significantly  less 
computer  storage  than  the  file  of  coefficients  produced.  Thus,  it  can  be  worthwhile 
to  generate  the  data  as  needed  rather  than  store  it  explicitly.  Such  generators  are 
easily  incorporated  in  GNET.  (This  approach  does  not  entirely  avoid  arc  length  infor- 
mation for  capacitated  problems  since  a record  must  be  maintained  for  each  arc  out  of 
the  basis  at  its  upper  bound.)  These  advantages  also  invite  development  and  use  of  a 
random  problem  generator  subroutine  to  replace  the  cumbersome  problem  fi ?es  and  cost 
of  using  NETGEN  [39]  to  generate  very  large  problems.  The  design  of  GNET  is  also 
consistent  with  the  use  of  explicit  arc  arrays  stored  in  peripheral  devices. 

Perhaps  the  most  important  potential  for  the  pure  network  solution  speed  of  GNET 
lies  in  more  general  large  scale  models  with  imbedded  networks  in  their  structure. 

The  multicommodity  distribution  system  design  code  of  Geoff rion  and  Graves  [22]  has 
been  revised  to  incorporate  GNET  to  repeatedly  solve  the  (thousands  of)  network  sub- 
problems. A goal  programming  code  for  network  problems  with  quadratic  objective 
functions  has  been  successfully  built  by  the  authors  with  GNET  used  as  the  key  sub- 
routine. 

The  theme  of  GNET  is  replacement  of  arithmetic  primal  simplex  computations  by 
simple  but  equivalent  logical  tests.  This  is  reminiscent  of  the  motives  for  generalized 
upper  bounding  and  suggests  that  network  factorization  may  prove  to  be  a viable  com- 
petitor for  GUB  in  general  linear  programming  systems  (31,  32). 

Lee  [43]  indicates  that  truly  huge  models  may  be  solved  by  mathematical  aggrega- 
tion and  develops  a wide  class  of  network  aggregation  methods  producing  surrogate 
problems  that  are  pure  networks  which  are  smaller  in  size  and  which  preserve  and 
exploit  special  global  structure  in  the  original  problem.  This  approach  is  intriguing 
due  to  the  curious  general  improvement  in  performance  encountered  when  solving  real 
models  rather  than  random  test  problems  of  equivalent  size. 
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